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Planet-disk interactions, where an embedded massive body interacts gravitationally with the 
protoplanetary disk it was formed in, can play an important role in reshaping both the disk and 
the orbit of the planet. Spiral density waves are launched into the disk by the planet, which, if 
they are strong enough, can lead to the formation of a gap. Both effects are observable with 
current instruments. The back-reaction of perturbations induced in the disk, both wave-like and 
non-wavelike, is a change in orbital elements of the planet. The efficiency of orbital migration 
is a long-standing problem in planet formation theory. We discuss recent progress in planet-disk 
interactions for different planet masses and disk parameters, in particular the level of turbulence, 
and progress in modeling observational signatures of embedded planets. 


1. Introduction 


The subject of planet-disk interactions has arguably en- 
tered its fifth decade (Goldreich and Tremaine|1980). The 
interplay between a massive perturber and a gaseous disk, 
and the resulting orbital evolution of embedded planets, 
nevertheless remains a rich subject in which significant 
progress has been made since the last Protostars and Planets 
chapter (PPVI, |Baruteau et al.12014). In addition, the sub- 
ject has been propelled forward through new observations 
of both fully formed planets and protoplanetary disks. Be- 
low, we first provide some relevant highlights of progress in 
both of these observational areas since PPVI. 


1.1. Planet population 


The population of known extrasolar planets has grown 
dramatically since the last Protostars and Planets review 


(Baruteau et al.|2014). We now have more than 4, 000 con- 
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firmed planets, amongst which around our nearest neigh- 
bour star (Anglada-Escudé et al.|2016), and this wealth of 
data reveals interesting patterns that may have bearing on 
planet formation and disk-planet interactions. Of particular 
interest are the radius valley (Fulton et al.|2017), and the ec- 
centricity dichotomy (e at POT chere ince lane 
systems are observed to have higher eccentricities than mul- 


tiplanet systems. These structures provide valuable links 
between the current architecture of planetary systems and 


their formation (e.g. 2017; (Ginzburg et al. 
2018| |Poon and Nelson|2020). A more direct link to for- 


mation can be made by observing exoplanets around very 
young, pre-main-sequence, stars. Such young planets have 
been detected by the radial velocity method around weak- 
line T-Tauri stars V830 Tau (Donati et al.]|2016), TAP- 
26 (Yu et al.|2017), and using transits for K2-33 
et al.|2016), V1298 Tau and AU Mic 
(Plavchan er al 2020) 
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Other new developments include characterizing exo- 
planet atmospheres, again with a possible link to their for- 
mation (e.g. 2018), the ability of GAIA to 
weigh very young planets (Snellen and Brown\2018), and 
the fact that the time baselines for exoplanet detections 
are now long enough to give information on planet oc- 
currence beyond the snow line (Fulton et al.[2021). This 
ever stronger link between observed planetary systems and 
their formation makes the study of disk-planet interactions, 
which will determine the orbital architecture of planetary 
systems before the protoplanetary disk disappears, ever 
more relevant. 


1.20. Disk observations 


Since PPVI, spatially resolved observations of proto- 
planetary disks have witnessed significant progress. De- 
tailed reviews of this progress can be found in the chapter 
led by M. Benisty, J. Bae, and also [Andrews] (2020). Here 
we briefly introduce disk observations that are most relevant 
to detecting signatures of disk-planet interactions. 

Today, high fidelity imaging observations of protoplan- 
etary disks on the global scale with the highest angular 
resolutions (S 0.1") and sensitivities are mainly carried 
out at two spectral windows. At optical to near-infrared 
(ONIR) wavelengths, a few 8-m class ground based tele- 
scopes equipped with adaptive optics (AO) and stellar halo 
suppression techniques have enabled us to image disks in 
both total and polarized intensity. A few representative in- 
struments of this category include the Spectro-Polarimetric 
High-Contrast Exoplanet Research (SPHERE) on the VLT 
(e.g., Benisty et al.|2015), the Gemini Planet Imager (GPI) 
on Gemini (e.g.,|Follette et al.|2017), and the Subaru Coro- 
nagraphic Extreme Adaptive Optics (SCExAO) on Subaru 
(e.g.. Currie et al.|2019). The diffraction limited angular 
resolution — ^40 mas, or < 6 AU at a typical distance of 
140 pc to the closest star forming regions — can usually 
be achieved in the H-band. At these wavelengths disks are 
usually extremely optically thick out to ~ 100 AU. Such 
observations probe scattered starlight from the disk surface 
determined by the spatial distribution of small dust parti- 
cles, ~ jym-sized or smaller. Such particles tend to be well- 
mixed with the gas and have an extended vertical distribu- 
tion instead of settling to the disk midplane. 

The second spectral window is at millimeter (mm) to 
centimeter (cm) wavelengths, which has been revolution- 
ized recently largely thanks to the Atacama Large Millime- 
ter/submillimeter Array (ALMA) and the upgraded Very 
Large Array (VLA). Dust continuum and gas line emission 
can be imaged by the two interferometers. With the former, 
we can trace the spatial distribution of dust particles of sub- 
mm to cm sizes, as they dominate thermal emission at these 
wavelengths. The latter type of observations allows us to 
probe the spatial distribution, temperature, and kinematics 
of gas in disks. Angular resolutions up to ~ 20 mas, or 
3 AU at a typical distance, can now be routinely achieved 
since the completion of the ALMA long baseline configu- 
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rations (ALMA Partnership et al.|2015). 
1.3. Disk properties 


Protoplanetary disks are complex entities, with differ- 
ent physical processes important at different locations (see 
chapter by Lesur et al.). Here, we start from very simple 
disk models, and summarise the main disk parameters that 
are important to understand disk-planet interactions. This 
sets the scene for discussions of recent progress using more 
realistic disk models in section[3] 

Protoplanetary disks consist mostly of gas, and the rel- 
evant equations describing their dynamics are therefore 
given by the equations of fluid dynamics: conservation of 
mass, momentum and energy. These need to be closed by 
an equation of state, expressing pressure in terms of den- 
sity and energy. If a barotropic equation of state is adopted, 
p = p(p), then there is no need to solve the energy equation. 
Such approximations are still widely used, in particular the 
locally isothermal approximation, where the temperature is 
a prescribed function of radius. The simple, locally isother- 
mal disk model can be made more realistic by including 
additional physics in the form of thermal effects (heating 
and cooling), magnetic fields (including non-ideal MHD) 
and solid particles in the form of dust. 

Disks are assumed to be close to vertical hydrostatic 
equilibrium, so that the vertical component of stellar gravity 
is balanced by a pressure gradient. If the disk is vertically 
isothermal, this leads to a Gaussian density profile with a 
pressure scale height H that is determined by the tempera- 
ture (or sound speed, cs): 

H = esf, (1) 
where X) is the local Keplerian angular velocity. Proto- 
planetary disks are usually thin, so that H/r = h < 1. 
This allows in some cases for a simplification through ver- 
tical averaging, after which we can effectively work in two 
dimensions and with a surface density instead of volume 
density. 

Disk masses are typically much smaller than the stel- 
lar mass, and the Toomre Q parameter is much larger than 
unity for large parts of the lifetime of the disk. This means 
that self-gravity can be safely neglected. This is not true 
for younger disks, which are thought to go through a phase 
where they are self-gravitating and potentially form giant 


planets by gravitational instabilities (e.g. 1978 


1997). It should also be noted that in some important 
cases self-gravity can have important effects when Q > 1 


but Qh ~ 1. This includes the Type I migration rate 


(Baruteau and Masset|2008b), global edge modes induced 
by planetary gaps (Lin and Papaloizou|2011), and vortices 
in low-viscosity disks (Zhu and Baruteau\2016} 
2021). 


A small percentage of the mass of the disk is in solid 
form, which is usually called dust. The canonical value for 


the dust-to-gas ratio is 1/100 (Draine|2011), which might 


indicate that dust is dynamically unimportant. However, 
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dust tends to cluster in pressure maxima, so that the local 
dust-to-gas ratio may be much higher than 1/100. This is 
where dust can in principle play a role in disk-planet inter- 
actions. 

Finally, angular momentum transport in disks is usually 
modelled by including a turbulent viscosity, usually through 


the a-prescription (Shakura and Sunyaev\1973): 


v = acs4H, 


(2) 


where v is the kinematic viscosity and a is a dimensionless 
constant. The underlying idea is that while molecular vis- 
cosity is many orders of magnitude too small to play a role, 
the disk is likely to be turbulent, with a maximum eddy size 
of ~ H (to fit into the disk) and maximum turbulent veloc- 
ities of ~ c, (otherwise strong damping by shocks). One 
then expects a < 1. While there are obvious limitations to 
modelling turbulence this way, the level of viscosity serves 
as a useful boundary for various regions of disk-planet in- 
teractions (see section[). 


1.4. Outline 


The rest of the chapter is organised as follows. We be- 
gin in section [2] by providing a basic overview of the vari- 
ous regimes of disk-planet interactions. We then go on to 
discuss areas of progress since the last Protostars and Plan- 
ets conference in section B] We discuss the tools used to 
make progress in section|4] In section|5]we discuss possi- 
ble observational signatures of disk-planet interactions, and 
conclude in section[6] 


Planet Mass 


Aerodynamic Coupling 


0 Disk Turbulent Viscosity 


Fig. 1.— Organization of basic planet migration regimes from 


McNally et al.\(2019a). 


2. Overview of Planet-Disk Interaction Regimes 


The interaction of a body orbiting a young star in the 
protoplanetary disk has a different character as a function of 
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the parameters of the solid body, the orbital trajectory, and 
the disk. To establish a rough orientation in the dominant 
physical regimes and the scales between them we use the 
map in Figure [I] 

We first summarize the basic flow structures in the disk 
resulting from an embedded planet. We then go on to de- 
scribe the various regions of Figure[I]that were well known 
at the time of the last Protostars and Planets review. 


2.1. Flow structure around an embedded object 
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Fig. 2.— Surface density perturbation around a Neptune-mass 
planet (q = 1074) on a circular orbit. Top panel: global overview 
of the disk, showing the prominent one-armed spirals. Bottom 
panel: zoom-in on the planet, with schematic streamlines (in the 
frame rotating with the planet) superimposed, indicating horse- 
shoe trajectories as well as trajectories responsible for the wave 
torque. 


The response of a gaseous disk to a massive embedded 
object can be divided up into a wavelike and a non-wavelike 
response. The wavelike response can be thought of to be 


generated at Lindblad resonances (Goldreich and Tremaine 
1979) and results in a one-armed spiral due to constructive 


interference of Fourier modes (Ogilvie and Lubow|2002). 


An example of a Neptune-sized planet (mass ratio q = 
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M,/M, = 10~*) embedded in a disk with H/r = 0.05 is 
shown in Figure 2| showing the prominent one-armed spi- 
ral density waves in the top panel. Additional spirals may 
emerge at larger distances from the planet, depending on the 
disk and planet properties (Bae et al.[2018). Such spiral pat- 
terns have been observed in a number of disks in scattered 
light as well as dust and gas emission (see section[5]. 
Close to the orbit of the planet, there is a wave-free zone 
if the disk is not self-gravitating and only weakly magne- 
tised (e.g. 
. This is where the 


fluid executes horseshoe turns in the frame corotating with 


the planet (e.g.|Dermott and Murray|1981). In a gas disk, 
the width of the horseshoe region depends both on the mass 
of the planet and the temperature of the disk (Paardekooper 


and Papaloizou|2009b). For low-mass planets, the half- 
width of the horseshoe region z,/ry ~ 4/q/h, while for 


high-mass planets the gas-free result is obtained, z,/r, ~ 
q'/3. The connection between the two mass regimes is such 
that for intermediate-mass planets the horseshoe region is 
wider than would be expected based on the low-mass case 
(Masser er al 2006) provide 
an expression of the width of the horseshoe region suitable 
in this transitional regime. In the bottom panel of FigureD]a 
close-up of the planet is shown, with schematic streamlines 
for both wavelike perturbations and horseshoe turns. 


2.3. Gas Drag / Type 0 


Bodies smaller than planetesimals of ~ 10 km drift due 
to gas drag (see chapter by Drazkowska et al.). In this case, 
gravitational interactions do not play a role. This regime is 
indicated at the bottom of Figure|1| While for larger bod- 
les, gravitational interactions usually dominate over aero- 
dynamics, in specific cases planetary bodies may undergo 
Type 0 migration (Adams et al. [2009 
[Becker et al.|2021). This kind of frictional interaction is 
also relevant for highly eccentric or inclined orbits 


2012). 
2.3. Viscous Type I 


Moving up in planet mass above the aerodynamic cou- 
pling transition in Figure |1| gravitational interaction with 
the disk becomes the dominant source of torques. This 
is dominated by the wave torque 


1980 1993) and the corotation torque or horse- 
shoe drag (Ward|1991). Towards larger values of turbulent 


viscosity, this is the realm of classic, viscous Type I migra- 
tion Tanaka et al.{2002). It is characterised by 
a corotation torque that is kept unsaturated by turbulent vis- 
cosity (Masse7]2001| 
. This 


means that while the reservoir of angular momentum that 
coorbital material can give to the planet is finite, viscosity is 
strong enough to replenish the angular momentum before it 
runs out. Towards lower values for the viscosity, coorbital 
angular momentum is not replenished fast enough, which 
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leads to a reduced, saturated corotation torque. Towards 
higher planet masses, torques on the disk become strong 
enough to modify the density structure of the disk, leading 
to the opening of a deep annular gap, which is the realm of 
Type II migration. 

A significant amount of effort has been made over the 
years to better understand Type I migration and to come up 
with parametrizations that can be used for example in N- 
body simulations and planet population synthesis models. 


Early examples include |Korycansky and Pollack| (1993), 


who tackled the two-dimensional isothermal case, and 
(2002), who provided torque prescriptions 
for three-dimensional, isothermal disks. These have the 
following form for the torque on the planet: 


(3) 


q? 402 

T=C 72 PT pps 

where q is the mass ratio planet/star, h = H /r is the aspect 
ratio of the disk, X is the surface density, ry the orbital 
radius and €), the orbital frequency. All spatially varying 
quantities are evaluated at the location of the planet, as in- 
dicated with a subscript p. The dimensionless constant C 
depends on the background density and temperature pro- 
files and is typically of order unity and negative, implying 
inward migration. It should be noted that unless the viscos- 
ity is strong enough to directly affect the horseshoe turns, 
the corotation torque is given by (nonlinear) horseshoe drag 
rather than the linear corotation torque (Paardekooper and 
When the assumption of an isothermal disk is relaxed, 
the corotation torque starts to play a very prominent role 


(Paardekooper and Mellema|2006a 


2008a||Paardekooper and Papaloizou\2008}\Kley and Crida 
2008). Based on the local gradient of entropy in the disk, 


it can even drive outward migration if the cooling of the 
disk is sufficiently inefficient for the disk to behave almost 
adiabatically (e.g. [Masset and Casoli|2009) 
et al.|2010). The constant C in equation (3) then depends on 
the local gradients of density and temperature, the level of 
turbulence, which governs the saturation of the corotation 
torque, and the cooling time of the disk. Parametrizations 
of the torque on the planet in this regime were presented in 


(2010) and |Paardekooper et al.| (2011), 
and subsequently by (2017) who pro- 


vided torque formulae adapted to three-dimensional disks. 


2.4. Viscous Type II 


When the planet becomes sufficiently massive, the spi- 
ral waves emanating from the planet deposit a significant 
amount of torque in the disk (typically via shocks), carving 
out a low-density annulus (i.e., a planetary gap) along the 
planetary orbit. The classical model of Type II migration as- 
sumed the density inside the gap to be low enough such that 
Type I torques play no role, and no gas flows crossing the 
gap. Under these assumptions, the planet motion would be 
locked to the same drift velocity as the background viscous 
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accretion of the disk (e.g. |Lin and Papaloizou|1986). In 


other words, the planet is locked inside its own gap, which 
slowly accretes onto the star. 

However, many hydrodynamical simulations have ob- 
served considerable gas flows crossing the planetary gaps 


since the end of the last century (e.g. 
1999, 
Masset and Snellgrove|2001). Gas was found to eas- 
ily drift through the gap region on horseshoe orbits. Recent 
simulations on Type II migration also found that a massive 


planet with a gap can migrate independently of the disk ac- 
cretion, owing to the gap-crossing flow (e.g. |Duffell et al. 
2014} Dürmann and Kley|2015||2017). A planet migrates 


together with the gap because the planet keeps digging a 
gap at its position. Based on these hydrodynamical simu- 
lations, a new simple model of Type II migration was pro- 
posed by [Kanagawa et al.| (2018b), to replace the classical 
model. This new migration model is dependent upon the 
accurate gap structure revealed by recent hydrodynamical 


simulations (e.g. |Duffell and MacFadyen|2013| 
2014 2015a). Such recent progress on the 


nature of gap structure and Type II migration will be dis- 
cussed at length in Section 3.4] 

There is some ambiguity in the semantic question of 
what “Type II" actually means. In some cases, “Type II 
migration” has been used to mean that the planet migrates 
in lockstep with the viscous drift rate of the disk (indeed, 
that is how the regime “Type II" was defined when it was 
first coined, see [Ward|1997). Under this definition, “Type 
II" migration would never happen, according to our present 
understanding. However, in this review we define "Type II" 
to simply mean the regime in which the planet opens a sig- 
nificant gap, enough to change the planet's migration rate. 

At what planet mass does this transition occur? The 
criterion proposed by [Crida et al.| had been widely 
used for many years, which is a combination of the thermal 
condition and the viscous condition. However, this crite- 
rion is no longer considered accurate, especially for low- 
viscosity disks. The gap-opening criterion has since been 
modified by recent hydrodynamical simulations (e.g. 
fell and MacFadyen|2013) and by traditional analytical con- 
siderations (e.g. 
[gawa et al.|2015a). The presently-understood gap open- 
ing condition is the same as the traditional viscous crite- 
rion (e.g. [Lin and Papaloizou]I986] Rajikov|2002) derived 
from the balance between the angular momentum flux of 
three-dimensional spiral waves excited by the planet (e.g. 
and the angular momentum flux of the 
viscous accretion disk. The classical thermal condition (i.e. 
Ruin Z H) is found to be unnecessary, as previously sug- 


gested by (2002). Thus, the gap opening condition 
is expressed as (Duffell and MacFadyen|2013| 
et al.|2015a) 


1/2 
q > 5n3/2 (=) = 5n5/291/2. (4) 
Tb p 


Planet-Disk Interactions 


This gap opening condition indicates that even a Neptune- 
mass planet will open a gap in a disk with a low-viscosity 
of a X 107^ and a disk aspect ratio of h = 0.05. The 
gap opening by such a low-mass planet was also directly 
shown by numerical simulations by[Duffell and MacFadyen] 
(2013), although this gap opening had been excluded by the 
classical thermal condition. 

Even when the gap opening condition is satisfied, the gap 
does not become empty immediately but the gap surface 


density decreases with a power-law of q7? (e. g.|Duffell and| 
[MacFadyen[2013). As the gap becomes sufficiently deep, 
the gap depth becomes a steeper function of q 
2014] [Duffell et al.|2020). The low-density gas inside the 


gap causes gap-crossing flows and accretion flows onto the 
planet even in a case of a Jupiter-mass planet. It also enables 
the planet to migrate independently of the disk accretion. 

For a viscous disk, the gap-opening criterion can be 
considered as a competition between timescales; i.e. can 
the planetary torques open the gap faster than viscosity 
can refill the gap? In the inviscid case, a different con- 
dition to equation should apply. In the absence of 
viscosity, gap opening is a competition between plane- 
tary torques and migration; i.e. can the planet open the 
gap on a timescale shorter than its migration timescale? 
This criterion is known as the “inertial limit" (Ward and| 
[Malik er al |2015), which 
is shown in Figure|l|as the Feedback Mass. In principle, 
other timescales could be compared with the gap-opening 
timescale, e.g. the lifetime of the disk, but the inertial limit 
is most likely the relevant criterion in the inviscid case. 


2.5. Type III 


Viscous Type I migration is driven by asymmetries be- 
tween the inner and outer spiral wave (the wave torque) and 
asymmetries between the two horseshoe legs due to radial 
gradients in density and temperature. Torque prescriptions 
are usually derived by keeping the planet on a fixed orbit 
and measuring the steady state torque on the planet 
[Paardekooper et a[]]OTD. If the 
planet is migrating, this potentially creates another source 
of asymmetry between the two wakes and in particular the 
two horseshoe legs. 

The basic picture is illustrated schematically in Figure[3] 
While in the non-migrating case, both horseshoe legs would 
be symmetric, if the planet is migrating towards the left of 
the figure, all material not interacting with the planet will 
drift towards the right in the comoving frame. This leads to 
the structure depicted in Figure[3] On the bottom leg, unper- 
turbed disk material streams past the planet, executing only 
a single horseshoe turn. The top leg consists of coorbital 
material that migrates with the planet. 

This setup was first analysed in [Masset and Papaloizou| 
(2003). If the planet has partially depleted its coorbital re- 
gion, this will lead to a different torque from both legs even 
in the absence of global gradients. Moreover, when migra- 
tion is slow compared to the horseshoe libration time, the 
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Fig. 3.— Schematic flow structure of the coorbital region in the 
comoving frame of a migrating planet. Orbital radius varies hor- 
izontally, while the direction of orbital motion is vertical. The 
planet is migrating towards the left of the figure. 


resulting torque is proportional to the migration rate, open- 
ing up the possibility of a runaway process 
Artymowicz|2004). This was studied further 
in (2008alb|c). For a more in-depth discus- 
sion see 2007). 


3. Areas of Recent Progress 


In this section, we review areas where there has been 
significant progress and interest, particularly since PPVI 
(Baruteau et al.\2014). 


3.1. 


Most effects on the planetary torque arising from ther- 
mal diffusion within the gaseous disk, and those due to the 
release of energy into the surrounding gas by an accreting, 
low-mass planet, were not taken into consideration in the 
work reported in PPVI. Since then, it has been discovered 
that the diffusion of heat in the vicin- 
ity of a low-mass planet can have a large impact on the 
disk torque. The gas that tends to accumulate in the vicin- 
ity of the planet is compressionally heated. The heat thus 
produced tends to diffuse outwards, and the planet's neigh- 
bourhood is colder (hence denser, to maintain the pressure 
balance) than it would be if the gas behaved adiabatically. 
The difference in gas density with respect to the adiabatic 
case takes therefore the form of a dense peak centred on the 
planet. At small distance from the latter, diffusion occurs 
on a timescale much shorter than that of advection, and this 
peak has a central symmetry with a density that decays as 
rr being the distance to the planet. At further distances, 
the diffusion timescale increases and becomes larger than 
the shear timescale. Advection then takes over and distorts 
the perturbed region: for a planet on a circular, non-inclined 
orbit, two lobes appear downstream of the planet (therefore 
one behind the planet in the outer disk and one leading the 
planet in the inner disk). As the disk is generally slightly 
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sub-Keplerian, the planet is outside its corotation, so that 
the outer lobe is more pronounced and dominates the force 
exerted on the planet. Since this dense lobe is trailing the 
planet, its contribution to the planetary torque is negative: 
in a disk with thermal diffusion, inward migration is faster 
than it would be in a similar, adiabatic disk. Owing to the 
shape of the density difference with respect to the adiabatic 
case, [Lega et al.| (2014) have called this effect the cold fin- 
ger effect. 'The new torque contribution arising from the 
lobes has been called thermal torque (Masset\2017). In 
the case studied by (2014), where the planet 


does not release energy into the ambient disk, it is more 
specifically called the cold thermal torque. The typical size 
À of the lobes can be obtained by equating the diffusion 
timescale A? /x (where x is the thermal diffusivity) and the 
shear timescale Q~!. One obtains: 


Ac yx/9. 


In planet-forming regions of protoplanetary disks, radiative 
diffusion is by far the main source of heat diffusion, and the 
thermal diffusivity takes the form: 


(5) 


e 16(y — 1)oT? 
3p? (R/u)s ’ 


where ^y is the adiabatic index, c is Stefan's constant, T' 
the midplane temperature, p the midplane density, 7 the 
ideal gas constant, jj the mean molar mass and « Rosse- 


land's mean opacity (Paardekooper et al.|2011). Eq. (6) 
has been checked by numerical experiments (Jiménez and| 


and found to yield values of the thermal dif- 
fusivity within 20 to 30 96 of that measured in runs involv- 
ing the radial thermal diffusion of hot, radially narrow re- 
gions. Evaluating A in the planet forming regions (within 
10 au) of protoplanetary disks shows that it is usually a tiny 
fraction (~ O(10~+)) of the pressure length scale H. Ob- 
taining converged results for thermal torques in numerical 
simulations is therefore quite challenging and requires very 
high resolutions (Chametla and Masset|2021), which ex- 
plains why thermal torques have been unnoticed for so long 
in studies of planet-disk interactions. Despite the minute 
size of the region involved in thermal torques, they can have 
a large impact on the net torque. In particular, for very 
low-mass planets (typically below one Earth mass), they 
can be larger than the differential Lindblad torque by an or- 
der of magnitude (Masser[2017). Note that so far there has 
not been any study of thermal torques in turbulent disks. 
Thermal torques are typically established over one orbital 
period (Chametla and Masset|2021). Should the turnover 
time of turbulent eddies with size ~ A be shorter than this 
timescale, one could expect a sizable impact on the value of 
thermal torques. This question requires further study. 

A faster inward migration when thermal diffusion is ac- 
counted for is not the end of the story. Low-mass planets 
may accrete pebbles or planetesimals. When that happens, 
they heat the surrounding gas and a process similar to that 
outlined above takes place, except for a change of sign: 


(6) 
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with respect to the non-accreting case, a hot, underdense 


region develops, with a two-lobe shape (Benitez-Llambay 
et al.[2015), as illustrated in Figs.[4|and[5] In the usual case 


of a sub-Keplerian disk, the outer lobe dominates the net 
contribution to the torque. This contribution, dubbed heat- 
ing torque, is positive. If the planet is sufficiently luminous, 
it can reverse the torque on a low-mass planet, driving it 
outwards. Using linear perturbation theory, one can obtain 
the expression of the thermal torque, when the distance of 
the planet to its corotation £p is much smaller than the size 


of the thermal lobes (Masset|2017): 


L 
(z = i) Teg ACT.) 


—lz 
Tthermal = ie 7e 
y A 


where L is the planet’s luminosity and Le the watershed 
luminosity at which the net thermal torque changes sign: 


i= 4nGMyxp. (8) 
^Y 

When x, becomes comparable to A, the thermal torque no 
longer scales with x, and tends to saturate to a value of 
order V2r(y — 1)/7(L/Le — 1)YyráQ2g?h.-3 
land Masset|2021). (2020) have used high- 
resolution, three-dimensional calculations in the shearing 
sheet to assess the validity of linear theory. They found an 
agreement at better than the 10 % level in the limit of low 
luminosity and low thermal diffusivity, and departures from 
linear theory at high luminosity or high thermal diffusivity. 

Thermal effects also have an impact on the time evolu- 
tion of eccentricity and inclination. These effects are de- 
scribed in section 

Thermal forces are the dominant contribution to the 
force exerted by the disk onto the planet up to a few Earth 
masses. As the mass increases past this value, their relative 
contribution to the net force decays slowly, and becomes 
negligible above 10 — 20 Earth masses. The decay of ther- 


mal forces has been studied in detail by 
(2020), in a different context: that of a point- 


like mass that travels through a uniform three-dimensional 
medium with thermal diffusivity x, and that releases energy 
at a constant rate L. This setup corresponds to that of dy- 
namical friction in a gas studied by [Ostriker| (1999), with 
the sole modifications of the inclusion of thermal diffusion 
in the gas, and of the energy released by the travelling mass. 
Owing to the existence of an axis of symmetry, this setup 
can be numerically tackled with two-dimensional meshes, 
and allows for a much more extensive study of the param- 
eter space, and a higher resolution, than the more complex 
situation of a planet embedded in a stratified gaseous disk 
with a sheared Keplerian flow. 
(2020) find that the thermal force exerted on the travelling 
mass is in good agreement with that predicted by linear the- 


ory (Masset and Velasco Romero|2017,| Velasco Romero and 
2019) for masses up to a few 


— Xs 
GY? 


M. (9) 


Planet-Disk Interactions 


and is O(Mc/M) times smaller at larger mass. The value 
of the critical mass M, depends on the location in the disk. 
It typically ranges from a fraction of an Earth mass just be- 
yond the snow line to several Earth masses at 10 au. Inside 
the snow line, by virtue of Eq. (6). the critical mass can be 
of several Earth masses, owing to the drop of opacity. Al- 
though no systematic study of the decay of thermal forces at 
larger planetary mass has been undertaken for planets em- 


bedded in protoplanetary disks, the findings of 
(2021) for a planet on a circular orbit are com- 


patible with a decay law similar to that found for dynamical 
friction. Similarly, the decay law for planets on eccentric 
or inclined orbits should closely match that found in studies 
of dynamical friction, since for eccentricities or inclinations 
larger than A/H, the thermal response reduces essentially to 
a hot or cold trail, depending on the value of the luminosity 
(Chrenko zr al. 2017). The results of [Eklund and Masset 
(2017), albeit obtained at low resolution, seem to confirm 
this expectation. 

An additional complication for planets with masses in 
excess of Me is that the width of their horseshoe region is 
comparable to or larger than the size À of the thermal dis- 
turbance (Chametla and Masset\2021). The thermal lobes 
are therefore distorted with respect to the shape they have 
at lower mass. The impact of this effect on the magni- 
tude of the thermal torque is yet to be assessed. A more 
extreme example of the consequences of a complex flow 


in the planet’s vicinity has been found by |Chrenko and 
(2019), who performed numerical simulations 


at high resolution with a non-constant opacity, and found 
that the heating torque is then strongly time variable, induc- 
ing episodes of inward as well as outward migration. This 
arises from the advection of heat in the planet’s immediate 
vicinity, on tightly wound streamlines. This process induces 
a highly time-variable configuration of the thermal lobes. 


3.2. A Generalization of Dynamical Corotation Torques 


The ‘Viscous Saturation’ vertical division in Figure[I]de- 
notes where viscous diffusion of angular momentum in the 
disk can maintain the background gradient of disk prop- 
erties across the corotation region. In the absence of vis- 
cous diffusion, phase mixing of material in the corotation 
region removes the background gradients that would oth- 
erwise drive the corotation torque. In this case, the con- 
trast and asymmetry of the disk material trapped on librat- 
ing orbits determines the corotation torque. This is partic- 
ularly applicable to low-turbulence wind-driven disks (e.g. 
Lesur et al.|2014), and has received much attention in recent 
years. For low-mass planets, this is the realm of dynamical 
corotation torques (see Figure D. 

The important flow pattern in the coorbital region is sim- 
ilar to that of Type III migration, see Figure[2| In contrast to 
Type III migration, however, we are now dealing with low- 
mass planets that do not open up a partial gap that drives 


Type III migration (Masset and Papaloizou|2003). How- 


ever, horseshoe drag for low-mass planets is, in an isother- 
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Fig. 4.— Relative perturbation of temperature arising from heat release by a low-mass planet, at the midplane of the disk. The central 
star is located at (x, y) = (0, 0). The left plot shows the case of a disk without radial pressure gradient. The planet is then located exactly 
at its corotation radius (materialized by the dotted line), and the two-lobe pattern is essentially symmetric. The right plot shows the more 
realistic case of a slightly sub-Keplerian disk. The planet is then located outside of its corotation radius, and the two-lobe pattern is no 
longer symmetric. Successive isocontours differ by a factor of 2, while the outermost contour corresponds to a value of óT'/T that is 
1/ 64°” of the peak value obtained in these calculations. These results have been obtained using pairs of three-dimensional numerical 
simulations with identical parameters, one in which the planet has a luminosity L and another one in which it is non-luminous. Here, the 
pressure length scale is H = 0.05a, the characteristic size A is ~ 0.14H (Eq.[) and the luminosity is L ~ 2.3L., where Le is given by 


Eq. (8). 


mal disk, driven by a vortensity contrast (Casoli and Masset 
2009), where vortensity is defined as the vertical component 


of the vorticity divided by the surface density: 

= V x v| z 
»» 
All we need for a nonzero corotation torque is a contrast 
between the vortensity in the coorbital region and the un- 
perturbed local disk. Then the two streamlines indicated in 
Figure 2|will carry a different value of the vortensity, yield- 
ing a net torque on the planet (Paardekooper|2014). The 
resulting vortensity distribution is illustrated in Figure [6] 
where the left panel shows a disk where viscosity is strong 
enough to maintain the original vortensity gradient across 
the corotation region, and the right panel shows the inviscid 
case, where a vortensity contrast has been built up through 
a combination of migration and radial disk flow 
ler al 20186). 

The criterion for a vortensity contrast between the coor- 
bital material and the rest of the disk to be maintained is 
that any diffusion should take place on longer timescales 
than the libration timescale (Paardekooper|2014). There- 
fore, the dividing line between the realm of dynamical coro- 
tation torques and viscous Type I migration is exactly the 
criterion for corotation torques to become saturated. If sat- 
uration happens, this means that dynamical torques enter 
the picture. 

One important difference between Type III migration 
and migration driven by dynamical torques is the way the 
vortensity or density contrast can be achieved. For high- 


(10) 


mass planets, one can rely on the planet clearing a par- 


tial gap (Masset and Papaloizou|2003), or the planet can 


be placed on a very steep density gradient (Pepliriski et al. 
2008alb|c). While the latter option is available for low-mass 


planets as well, migration itself is capable of generating a 
vortensity contrast. Since the planet migrates together with 
all coorbital material, in the absence of diffusion the vorten- 
sity associated with the coorbital material will be the same 
as the disk vortensity at the starting location of the planet. 
Therefore, if the planet is migrating in a disk with a ra- 
dial gradient in vortensity, a contrast will build up automat- 
ically. The resulting dynamical corotation torque can either 
enhance or slow down migration. 

If the migration direction of the planet (driven by the 
total torque) is the same as the direction in which the coro- 
tation torque alone would push the planet, the planet will 
be pushed into a runaway migration, much as in the case 
of Type III migration. If, on the other hand, the migration 
of the planet goes against the will of the corotation torque, 
migration will slowly grind to a halt (Paardekooper|2014). 
Adding some viscosity leads to residual migration, driven 
by how much angular momentum is able to diffuse into the 
coorbital region. 

The picture can be generalised to non-isothermal or ra- 
diative disks. In this case, the thermal memory of the li- 
bration island gives rise to new migration regimes at high 


optical depths (Pierens|2015| Pierens and Raymond|2016 


Another generalization comes in the form of a magnet- 


ically braked disk (McNally et al.|2017). The poorly ion- 
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Fig. 5.— Relative perturbation of the density field arising from 
heat release at the midplane of the disk corresponding to the right 
plot of FigureH] The isocontours have broadly the same shape as 
in that figure. The sign of the perturbation is the negative of that 
of the temperature: the lobes are under-dense. 


ized planet-forming regions may well be virtually free of 
magnetic turbulence, except for a narrow layer near the top 


from which a wind is launched 
[sel et al. [2015). If also the Hall effect is active, this creates 
significant laminar horizontal magnetic fields and associ- 
ated Maxwell stress 
[Béthune et al.[2017). This stress generates a radial flow in 


the disk, and while a radial flow is expected in viscous ac- 
cretion disks as well, in the laminar magnetic case there is 
no associated turbulence or diffusion. This means that even 
non-migrating planets would experience dynamical corota- 
tion torques, as the radial disk flow sets up the vortensity 
difference between the coorbital region and the rest of the 
disk. Depending on the ratio of the natural migration rate of 
the planet and the speed of the laminar disk flow, dynami- 
cal torques can either force the planet to migrate at the same 


speed as the disk flow, or cause a runaway (McNally et al. 


2017\|2018b). It is worth noting that an extension to higher 
masses has been shown to work (Kimmig et al.|2020). 


Above studies were mostly done in two spatial dimen- 
sions, working with vertically averaged quantities. While 
many of the features of horseshoe drag carry over in three 
dimensions 
[2017). there are important differences. One is that while 
migration appears to be insensitive to surface winds, in ra- 
diative disks an additional negative torque arises that is as- 


sociated with buoyancy waves (Zhu et al.|2012c| |McNally 
et al.12020a). 


3.3. Intermediate Mass Planets and No ‘Inertial Limit’ 


to Planet Migration 


Increasing the planet mass while keeping the turbulent 
viscosity small brings us through the "Feedback Mass’ 
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boundary in Figure This is where planets are mas- 
sive enough to significantly alter the density profile in their 
vicinity. 

The emerging new paradigm where protoplanetary disks 
are largely laminar with almost no turbulent diffusion 
has sparked renewed 


interest in the concept of the "inertial limit’ of planet migra- 
tion (Ward and Hourigan|1989 
2002). In 1D models of laminar disks, all plan- 
ets will open up deep gaps unless they can migrate through 


the disk fast enough (see also 2015 
[Bitsch[2017). For increasing planet mass, the density per- 


turbation induced by the planet provides a negative feed- 
back and migration stalls at the ’inertial limit’ 
[Hourigan|1989). The root cause of this stalling is the asym- 
metry in the induced density perturbation in the inner and 
outer disks. 

Early explorations using multidimensional hydrody- 
namic simulations did indeed find that migration would 


stall or slow down near the inertial limit (Li et al.|2009 
2010). |Fung and Chiang (2017) found that super- 


Earths in inviscid disks could stall their migration, while 
still driving accretion onto the central star through the ac- 
tion of their density waves (Goodman and Rafikov|2001). 
Migration stalling in combination with pebble accretion 
can potentially explain the dichotomy between inner super- 
Earths and outer Gas Giants (Fung and Lee|2018). 

disks with low levels of viscosity come with their own 
set of challenges. It has been known that intermediate-mass 
planets can excite vortices at the edges of their partial gap 
(Koller er al 003). Y was found in|McNally etal. (2019) 
that such vortices can seriously disrupt the picture of the in- 
ertial limit. An example is shown in Figure[7| At the lowest 
resolution, r2, which has 23 zones per scale height, migra- 
tion halts at around 3000 orbits, which has been seen before 
in e.g. (2010). However, soon after this, vortices 
emerge around the orbit of the planet, and the interaction 
between the planet and the vortices pushes the planet back 
on an inward trajectory. Eventually, at around 7000 orbits, 
it enters a phase of rapid, inward Type III migration. While 
the average migration rate is slow compared to viscous Type 
I (indicated by the dotted curve), migration is inward with- 
out any prolonged episodes of stopping. Moreover, dou- 
bling the resolution (r4) gives different results, with an early 
Type III phase that makes the average migration rate only 
about a factor of 2 slower than viscous Type I. Doubling the 
resolution again leads again to a different trajectory, now 
without any Type III episodes. This indicates that numer- 
ical convergence in these types of inviscid simulations is 
hard to achieve (McNally et al.|2019a]. 

Towards larger values of the viscosity, we cross the ’vis- 
cous smoothing’ border in Figure [1] towards the realm of 
smooth feedback. This is where viscosity can partly smooth 
out the action of vortices, while feedback effects are still 


important. This typically happens at v 1077 rot (cor- 
responding to a ~ 10~° for typical values of H /r) 
2010! McNally et al.|2019a). Migration in this regime can 
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Fig. 6.— Comparison of the coorbital region for two planets at the same orbital location. Left panel: Viscously driven inflow case, 
the vortensity gradient across the corotation region is maintained close to the background disk value by the viscosity, resulting in an 
unsaturated corotation torque. Right panel: Laminar magnetic stress driven inflow case, where the vortensity of the librating streamlines 
evolves due to a combination of the planet radial motion and the history of magnetic torques acting on the librating material. Figure 


from|McNally et al.\(2018b). 


be very slow compared to viscous Type I migration. 


3.4. The Nature of Type II Migration in Viscous Disks 


The classical Type II migration rate assumes an empty 
planetary gap with no crossing flow. With an idealized deep 
gap, the planet is locked between the outer and the inner 
disks and must migrate together with the radial viscous evo- 
lution of the disk. However, since around the turn of the 
century, many hydro-dynamical simulations had found sig- 
nificant gap-crossing flows for deep gaps formed by Jupiter- 
mass planets, as mentioned in Section After PPVI, 
highly accurate hydrodynamical simulations have revealed 
structure of the planetary gaps and Type II migration rates 
for very wide parameter ranges of the planet mass Mp, the 
disk viscosity v, and the disk scale height H. Owing to the 
reliable numerical results, a new qualitative and quantitative 
picture of the gap structure and Type II migration have been 
constructed. This recent progress is reviewed below. 


3.4.1. Gap Depth and Gap Profile 


Structures of almost axisymmetric gaps are described by 
radial profiles of the azimuthally averaged surface densities 
of disks. One of the most important quantities characteriz- 
ing a radial profile is the minimum surface density in the 
gap, Ueap, or the ratio of X4, to the unperturbed value, 
Xo, which is called the gap depth. Usually, the radial pro- 
file attains its minimum surface density at the planetary or- 
bital radius, ry if the high-density region around the planet 
is excluded at the density evaluation. Another important 
quantity for comparison to observations is the gap width. 

An intensive survey of accurate hydrodynamical simu- 
lations on the gap structure was first carried out by [Duffel] 


and MacFadyen| (2013) for wide parameter ranges of Mp 


and v. Their simulations with analytical considerations in- 


10 


dicate that the gap depth is given by 


Yap 


= 29)-! 
xe -(rK/9) ^, 


(11) 


where K = q?/(h?o) is a non-dimensional parameter 


(Kanagawa et al.|\2015a). This empirical formula for 


the gap depth is consistent with the numerical results by 
(2004). The formula was also confirmed by 
later surveys of hydrodynamical simulations and explained 
by zero- and one-dimensional analytic models (e.g. 
2017}. The gap open- 


ing condition (4) is consistent with equation (11) because it 
is rewritten as 3,45 /X9 < 0.5, using E 

Some studies proposed more accurate models of the gap 
depth than equation (11) for deep gaps of K > 10? (e.g. 
Duffell[2020). Y 
is useful to introduce another gap surface density, X54, av, 
which is radially averaged over the radial range between 
r+ = r + 2Ryin. The radially averaged surface density is 
physically more meaningful than X44, (c U(rp)), because 
the planet interacts mainly with the gas at r+ rather than 
at ry. Equation has been found to be valid even for 
K > 10°, if it is used for Xsap,av instead of Ueap. Accurate 
empirical models for the gap width and the radial profile of 


X(r) have also been constructed (e.g. [Duffell|2015a] 2020] 
[2017), including public tools avail- 
able for download to compute the complete density profile 
X(r) for given planet and disk parameters 
see Figure[8). 

For very massive planets with q > 3 x 10^? or for disks 
with a low viscosity, structures of the deep gaps are unsta- 
ble and the disks become eccentric (e.g. 
2014). The above gap models are not applicable to 
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Fig. 7.— Migration of a q = 1.25 x 10^? planet in an inviscid 
disk with H/r = 0.035. The three resolutions have 23 (r2), 47 
(r4) and 93 (r8) grid cells per scale height at r = rp. Figure from 


McNally et al.|(2019a) 


such eccentric disks. 


3.4.2. Type II Migration Rate in Viscous Disks 


The classical Type II migration rate is given by (e.g. [Lin] 
and Papaloizou\1986\|Syer and Clarke|1995) 


3v 2 
Tar (Mp < AXors); 
Ur,classical = 2 
; 3v AXor5 
—— M, > AXor2). 
2v M, Voa AX orp) 


(12) 


In this equation, the former is called the disk-dominated 
case and the later is the planet-dominated case. The con- 
stant A is chosen to be 2-10, depending on the authors (e.g. 
Baruteau et al.\2014). 

To verify the classical migration model, many hydro- 
dynamical simulations measured Type II migration rates, 
varying the disk parameters and the planet mass (e.g. 


2014; |Dürmann and Kley|2015 
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Fig. 8.— Comparison between the gap model of (2020) 
and a wide range of numerically calculated gaps. Black dashed 
curves are the analytical model, and colored solid curves are the 
numerically calculated gap profiles in steady-state disks. The fidu- 
cial disk model used q = 5.1 x 1074, a = 0.01, h = 0.05. Each 
panel picks one of these parameters and varies it, keeping the oth- 
ers fixed at their fiducial values. Figure from[Duffell] (2020) 


2018b 2018). Those results suggest that the 


migration rate can vary as v, c VXo/M, for relatively 
massive planets (although an anomalous additional depen- 
dence on ^ may still be present, e.g. [Duffell et al.|2014). 
This dependence is reasonably consistent with the planet- 
dominated case of the classical model. 

However, several points were found to be starkly incon- 
sistent with the classical model. For a relatively massive 
disk, a gap-forming planet can migrate faster than the radial 
flow velocity of the viscous accretion disk, 3v /(2rp). The 
observed migration rate also depends on the disk aspect ra- 
tio h for a fixed v. These results cannot be explained by the 
classical model. showed that no pile- 
up of the outer disk or no depletion of the inner disk occurs. 
This means that a considerable gap-crossing flow exists. 


A new picture of Type II migration proposed by [Kana-] 
(2018b) appears to be much more consistent 


with numerical results. In this picture, the gap-forming 
planet migrates with the Type I torque determined by the 
surface density at the gap bottom instead of the unperturbed 
one. Then, using equations and (11). the migration rate 
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where the gap opening condition (i.e. a very large K) 
is assumed and C is set to be ~ —3 at the second equality. 
The new migration rate resembles that of the classi- 
cal migration rate in the planet-dominated case, except the 
dependency on the aspect ratio h. This new model of the 
migration rate can also describe the continuous transition 
from the Type I migration to Type II. 

The new migration model has been used in population 


synthesis calculations for giant planets (e.g. [Ida et al.|2018 
2019 2021). Note that the new 


Type II migration rate is further reduced by a rapid gas ac- 
cretion onto the planet (Tanigawa and Tanaka\2016} 
Tanaka et al. [2020] [Wang et al.]2020). When 
rapid gas accretion onto the planet is regulated by a rela- 
tively low disk accretion rate, the planetary accretion further 
reduces the surface density in the gap, Xap. Thus the short 
supply reduces equitably the planetary migration rate and 
accretion rate since both are proportional to Xa 
[et al. |2020). This additional slow-down of Type II migra- 


tion should also be included in population synthesis calcu- 
lations. 

Further attempts have been made to refine the new Type 
II migration model. (2019) demonstrated 
that the back reaction of the net torque causes imbalance in 
the surface density between the outer and inner disk, which 
alters the torque imbalance (i.e. the net torque) on the planet 
in turn. They showed that this positive feedback enhances 
the net torque particularly in the low-viscosity case. The 
torque imbalance due to the outer gap edge heating by the 
central star was examined by two- and three-dimensional 
simulations (Hallam and Paardekooper[2018 
Nesvorny|2020). Further debate on the new Type II migra- 
tion model is ongoing 
Tanakal2020). 


One element still missing from existing Type II migra- 
tion models is the allowance for outward migration, which 
can occur due to a number of mechanisms. Migration (and 
growth) of super massive planets with q > 107? (transition- 
ing from a planet to a binary) was investigated by Duffell et 
al. (2020b). They found outward migration for g > 0.05 
as well as significant accretion onto the secondary, always 
driving the system toward an equal-mass binary. Similar 
outward migration of supermassive planets was reported by 
[Dempsey et al.|(2021). They posited that the outward mi- 
gration is due to the unstable eccentric outer disk. Interest- 
ingly, outward migration is not predicted by classical Type 
II migration, nor is it predicted by the new migration formu- 


las o (2018b). This by itself suggests there 


is substantial room for improvement of analytical models. 
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For intermediate planet masses, other effects can con- 
spire to yield, at least temporarily, outward migration. In 
disks with very low viscosity, super-Earths of q ~ 1075 
can open up gaps that are unstable to the Rossby Wave 


Instability (Lovelace et al.||1999). The interaction of the 


planet with the resulting vortices can drive outward migra- 
tion (McNally et al.|2019a), similar to what was seen for 
Saturn-mass planets Tt paseo Outward 
migration of super-Earths was also seen by (2015b) 
in disks of moderate viscosity. 


3.5. 


Most work so far has concerned the gas component of 
the disk only. While the dust component is recognized to be 
important for observations (e.g.|Paardekooper and Mellema] 
see also section 5), its effect on disk planet interac- 
tions has largely been ignored due to the small dust mass 
compared to the gas. Recent work, however, has recog- 
nized that the solid component can produce two new torque 


effects on the planets, the dust scattering torque 
Llambay and Pessah|2018}|\Chen and Lin|2018) and aero- 
resonant migration (Storch and Batygin|2019). 


The dust scattering torque arises from dust particles that 
stream past the planet. The resulting torque bears resem- 
blance to Type III migration, where the dust particles fol- 
low streamlines as depicted in Figure[3] The strength of the 
Type III torque is such that it can often compensate for the 


small dust mass compared to the gas (Benitez-Llambay and 
2018). This effect is most prominent for larger parti- 


cles that are not very well coupled to the gas. |Chen and Lin 
(2018) considered smaller particles, and found oscillatory 


torques from the coorbital region for high dust-to-gas ra- 
tios, together with vortex formation. In the thermodynamic 
framework of (2017), the vortices can be 
attributed to baroclinic effects stemming from a variation 
in dust-to-gas ratio (Chen and Lin|2018). More recently, 
(2020) found that a dusty version of dynam- 


ical corotation torques (see section can significantly 
slow down migration. 

Aero-resonant migration in- 
volves solid bodies of sizes between 10 m and 10 km. When 
these drift inward towards a planet due to gas drag, they 
get caught in mean motion resonances. Here, they can ex- 
change angular momentum with the planet and drive inward 
migration. The resulting migration rate can be comparable 
or larger than the viscous Type I migration rate, in particu- 
lar in the inner disk (Storch and Batygin|2019), if the total 
mass of planetesimals captured in mean motion resonance 
consists of a significant fraction of the mass of the planet. 

Gap-opening planets were considered by 
(2019), who found that dust feedback on the gas decreased 
the gas surface density at the outer edge of the gap, thereby 
reducing the torque on the planet. For wide and deep gaps, 
the direction of migration can reverse due to dust feedback 


(Kanagawa[2019). 


Impact of Embedded Solids on Planet Migration 
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3.6. Migration in Wind-driven Protoplanetary Disks 


The realization that protoplanetary disks could be largely 
laminar, with accretion driven by a wind rather than turbu- 
lent diffusion has 
led to new studies of disk-planet interactions in such disks. 
Consequences of the low level of turbulence have been dis- 
cussed already in sections [3.2] and [3.3] Here, we focus on 
recent progress specific for wind-driven disks. 

At low planet masses, the surface accretion flow driven 
by a wind is likely largely decoupled from the planet mi- 
gration (McNally et al.[2020a). Any surface accretion flows 
are too high up in the disk, where gas densities are too low, 
to influence the migration speed of the planet. No verti- 
cal communication between the accreting layers and the 
midplane was observed in 3D simulations 
P0203). 

Most studies of disk-planet interactions in wind-driven 
disks use a vertically integrated approach, where the effect 
of the disk wind is incorporated as an additional torque on 
the disk. incorporated the model 
of (2010) into a one-dimensional viscous evo- 
lution model for the disk, and constructed migration maps 
based on the formulae of (2011). Type 
I migration was found to be strongly reduced for a wide 
range of wind parameters due to the modified surface den- 
sity slope induced by the disk wind. Effects of rapid, wind- 
driven gas flow on the saturation of the corotation torque 
was investigated in [Ogihara et al.|(2017). Systems of low- 
mass planets in wind-driven disks were found to undergo 
dynamical instabilities, thereby avoiding resonant configu- 

At higher planet masses, 2D simulations of higher-mass 
planets in wind-driven disks showed that dynamical torques 
play a role in this mass regime (Saturn-Jupiter) as well 


(Kimmig et al.|2020). Fast enough accretion flows can drive 


planets into Type III migration that can be directed outward 


(Kimmig et ai.|2020). 


3.7. Migration of Planets Formed by Gravitational In- 
stability 


Gravitational Instabilities (GIs) may form planets from 
collapsing clumps in very early massive disks 
[Boss||1997). The physical regime 
is quite different from planet migration in a fiducial GI sta- 
ble disk. It was studied in detail in|Baruteau et al.|(2011), 
who found that newly formed fragments would migrate in 
the Type I regime, despite being Jupiter-sized. This is be- 
cause the time scale of GI to form a clump is the orbital 
time scale, which is much shorter than the gap formation 
time scale. The planet therefore has no time to open up a 
gap and starts Type I migration until it slows down enough 


to open a gap (Baruteau et al.|2011] |Malik et al.|2015). 


Interaction with global spiral modes can slow down or even 


halt migration (Michael et al.|2011). 


Studies with more realistic thermodynamics in the form 


of radiative disks were presented in (20122), |Sta-| 
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(2015), |Vorobyov and Elbakyan| (2018) and 
matellos and Inutsuka| (2018). They found that after an 


initial phase of fast migration, planets are able to open up 
gaps and slow down or even reverse the direction of migra- 
tion. The radiative feedback of accretion onto the planet 
is important in determining whether the planet can survive 


at large orbital distances (Stamatellos|2015). Recently, the 


importance of thermodynamics was further emphasized in 
[Rowther and Meru} (2020), who found that migration halts 
if the inner parts of the disk are gravitationally stable. Sub- 
sequent tidal and thermal evolution of clumps can change 
their mass, size and dust content, such that the resulting 
planet population could show a great deal of diversity (see 
e.g. [Nayakshin 2017). 

A code benchmarking study on planet migration and 
growth in self-gravitating disks was presented in [Fletcher] 
(2019). Migration rates between the various codes 


varied between 10% and 50%, depending on the numerical 
setup. In particular, they pointed out the importance of the 
artificial viscosity prescription and the way accretion was 
handled in the simulation in determining the migration out- 
come. 


3.8. Migration near the inner disk edge 


Here we discuss recent progress on disk-planet interac- 
tions near the inner edge of the protoplanetary disk. If the 
disk is truncated by the stellar magnetic field, there will be 
a low density inner magnetospheric cavity connected to the 
bulk of the disk. The steep, positive density gradient at the 
edge of the cavity can have important effects on migration. 

If the width of the cavity edge is smaller than the width 
of the horseshoe region, an approximate expression for the 
corotation torque is the one sided horseshoe drag 
2017). The resulting outward migration as the cavity ex- 
pands is able to break up resonances formed by earlier 
convergent migration ( magnetospheric rebound’, [Liu et al.] 
2017) [Liu and OrmeIQOT7) 

The question if and where planets are stopped near the 
inner edge has received a lot of attention. Wave reflec- 
tion off sharp edges was known to affect the Type I torque 
(Tsang|2011), and it was shown in|Miranda and Lai|(2018) 
that this can lead to the planet stopping as far out as 3 times 
the radius of the inner edge. It was argued in [Romanova] 
that disks may have several edges where plan- 


ets could stop migrating. The migration of a resonant chain 


of planets into the cavity was studied in |Ataiee and Kley 
(2021), who found that even chains of planets are efficiently 
stopped, and found that overstable librations (see|Goldreich 


and Schlichting|2014) may give rise to compact systems. 
Flock et al. (2019) studied the effect of the silicate sublima- 


tion front on the location where planets halt their migration, 
and found that Earth and super-Earth planets stop at periods 
ranging from 10 to 22 days, in good agreement with obser- 
vations. 

Near the inner disk edge is also the place where star- 
planet tidal interactions might play a major role in driving 
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orbital evolution. In particular, planets that stall their mi- 
gration near the inner disk edge may undergo further tidal 
migration once the disk has disappeared. For an overview of 
star-planet tidal interactions we refer the reader to 


(2014 2020). 


3.9. Eccentricity and Inclination Evolution 


Planet-disk interactions alter the planet’s orbital eccen- 
tricity and inclination as well as its semimajor axis, and 
their evolution should be treated consistently 
[2020]. For planets that are unable to open a gap in the disk, 
small eccentricities and inclinations e,i < h lead to small 
horizontal and vertical excursions within the disk. Until re- 
cently, it was believed that such e and ? would be damped on 
a timescale much shorter than that of Type-I migration, ow- 
ing to the emission of density and bending waves in the sur- 
rounding disk (Tanaka and Ward|2004). However, different 
conclusions have been reached by recent numerical simula- 


tions (Eklund and Masset|2017| |Chrenko et al.|2017) and 
linear analysis (Fromenteau and Masser[2019) that include 


radiative thermal diffusion and treat the accreting planetary 
embryo as a heat source that produces a hot, underdense 
trail in the disk (see Section B1). as the associated forces 
are typically larger than those due to the emission of waves 
for planets of mass comparable to or lower than the criti- 
cal mass of Eq. (9). If the planet's luminosity is below the 
critical value of Eq. (8). then e and ? are damped (and typ- 
ically much faster than in an isothermal or adiabatic disk); 
above the critical luminosity, both e and 7 grow to values 
comparable to h. The growth of e is significantly faster and 
tends to dominate the evolution, and the effect appears to be 
suppressed when the planet exceeds a few Earth masses. 

It has been appreciated for some time that there is a qual- 
itative change in behaviour when e = h, as this leads to 
supersonic relative motion between the planet and the disk. 
have proposed a simple prescription, based 
on dynamical friction calculations, for the damping of a, e 
and i of low-mass planets, that smoothly spans the sub- and 
supersonic regimes. 

For planets of about one Jupiter mass or above, capa- 
ble of opening a gap in the disk, the evolution of e and 
i over wide ranges of those quantities has been measured 
in numerical simulations by [Bitsch et al.] (2013), omitting 
the thermal effects discussed above. Generally e and i are 
both damped, although the damping timescales are greatly 
increased when these quantities are larger and the planet 
moves supersonically with respect to the disk. The mea- 
sured damping rates have been compared with theoretical 
expectations and used to explain observed systems 
ladis et al.|2017). Important exceptions to this behaviour 
are found (i) when e and 2 are both small and the planet 
is sufficiently massive to open a deep gap, in which case e 
grows; (ii) when 7 = 40°, in which case there can be non- 
monotonic behaviour related to Lidov-Kozai cycles. 

Earlier work on the eccentricity evolution of planets that 


open deep gaps (Goldreich and Sari|2003) had suggested 
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that e grows as a result of eccentric Lindblad resonances 
when e exceeds a small critical value, such that the oppos- 
ing eccentric corotation resonances are saturated 
land Lubow) 2003). measured 
€ > 0 in 2D numerical simulations for an intermediate 
range of e (0.04 < e < 0.07 for q = 107°), attribut- 
ing the upper limit to the collision of the planet with the 
disk. Eccentricity growth requires a clean gap to eliminate 
the coorbital Lindblad resonances, and is possible when the 
dimensionless parameter K = q?/ah? exceeds a thresh- 
old of 103—104. simulated cases with 
q = 1.3 x 107? and with two disk:planet mass ratios over 
300000 orbits, finding that both planet and disk become ec- 
centric. Differently from [Duffell and Chiang] (2015), they 
found that no seed eccentricity is required, and that e can 
grow to a value significantly larger than h. In 2D simula- 
tions with q = 10^? and very low viscosity, 
found excitation of e up to a much higher value of 
0.25, but the same did not occur in 3D simulations. 

Since gap-opening planets often have comparable angu- 
lar momenta to the disks in which they orbit, and since e 
and i are exchanged between the planet and the disk through 
secular gravitational interactions, it can be important to ac- 
count for the coupled motions of the planets and disks, 
taking into account the mean-motion resonances that con- 
tribute to excitation and damping of e and i. |Teyssandier| 
formulated a linear theory of this inter- 
action for small e, highlighting differences in the behaviour 
of e in 2D and 3D disk models. This type of analysis can 
explain the growth of disk eccentricity for planets above a 


few Jupiter masses (Teyssandier and Ogilvie[2017), and\Ra-| 
(2018) explained the long-term evolution of the 


eccentricities of the planet and disk in their simulations in 
terms of coupled modes of the system. 

(2019) proposed that the conserved 
angular-momentum deficit (AMD) of a system of two giant 
planets could be transferred from the outer planet to the in- 
ner one, conferring it with a large eccentricity, through the 
dispersal of the protoplanetary disk. However, the effect 
may be suppressed if the disk is also allowed to acquire 
(and to damp) eccentricity (Teyssandier and Lai|2019). 

If a planet migrates into a central cavity in the disk, then 
both the disk and the planet can become significantly ec- 
centric as a result of the eccentric Lindblad resonances that 


remain in the massive part of the disk (Teyssandier and 
2016). (2021) found that e as large 


as 0.4 could be achieved in this way. 

Several recent studies have addressed the inclination 
evolution of planets on circular orbits. |Arzamasskiy et al. 
measured the effect of non-zero 7 on the inclination 
damping rate and migration rate of low-mass (q = 1074) 
planets, finding that both increase for small non-zero i, but 


are greatly reduced for large i. (2018) and 


(2019) showed further that planets that open a deep 
gap break the disk into independently precessing inner and 


outer disks, resulting in a warped shape. 
In conclusion, it appears that there are several situations 


Paardekooper, Dong, Duffell, Fung, Masset, Ogilvie, Tanaka 


in which planet-disk interactions can cause a planet to ac- 
quire a non-zero orbital eccentricity. The inclusion of non- 
trivial thermal physics makes a significant difference to the 
outcome for low-mass planets. Since the eccentricity dy- 
namics of gap-opening planets can be strongly coupled to 
that of the disk, it may also be important in that context to 
include the effects of heating and cooling, along with 3D 
dynamics. 


3.10. Multiple Planets 


This section reviews work on planet-disk interactions 
in the generalized case of multiple planets, restricted to 
work where the details of the planet-disk interaction is high- 
lighted, including the formation of resonant chains. In 
particular the existence of resonant chains is an important 
observational diagnostic for planet-disk interactions (see 


also |Baruteau et al.|/2014). Particularly well-known sys- 


tems involving chains of mean motion resonances include 
Trappist-1 and Kepler-223 
(2016). In order to account for the diversity of orbital pe- 
riod ratios between adjacent pairs of planets, for example 
amongst the compact multi-planet systems discovered by 
Kepler, several studies in recent years have focused on how 
to break any resonant chains that arise during planet-disk 
interactions. 

Dynamical instabilities are a prime candidate for taking 
planetary systems away from resonance. This idea was first 


presented in (2012), and subsequent work 
(e.g. 2014 2017) showed that 


the distribution of planets as produced by disk migration 
followed by dynamical instabilities can match the observa- 


tions (/zidoro et al.|2021). 
Ogihara et al.|(2018| see also section [3.6) studied the 


formation of close-in super-Earths in wind-driven disks. In 
cases where the surface density profile is such that Type I 
migration is significantly reduced, it was found that dynam- 
ical instabilities set in, destroying any resonant chains that 
Type I migration had set up. The resulting distribution of 
period and mass ratios can reproduce the observed distribu- 


tion of close-in super-Earths (Ogihara et al.|2018). 


For a pair of planets in both viscous and inviscid Type 
I migration, (2018) found that while 
convergent migration creates mean motion resonances be- 
tween the planets, these can be destroyed by overstable li- 
brations, as suggested by|Goldreich and Schlichting (2014). 
This effect can potentially explain the relative lack of mean 
motion resonances in Kepler systems 
ider[20T8). 

The role of viscosity in forming resonant chains was ex- 
plored in (2019b). It was found that vis- 
cous disks indeed tend to form resonant chains, but that the 
picture changes qualitatively for inviscid disks. When the 
emergence of vortices (see section B.3) starts to make mi- 
gration paths non-smooth, convergent migration into reso- 
nances is disrupted. Overall, migration is still directed in- 
ward, creating compact planetary systems that are not pro- 
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tected by resonances. This frequently leads to dynamical in- 
stabilities in addition to closely packed, stable non-resonant 
systems (McNally etal 2019). 

Above studies indicate that while Type I migration al- 
most inevitably leads to resonant chains of planets, this is 
by no means the end of the story. Both viscous and inviscid 
disks have plausible mechanisms for breaking the chains, 
either before or after the gas disk has disappeared. 

The role of the disk mass was explored both in 
and[Ataiee and Kley|(2020).(Marzari|(2018) studied 
the effect of a finite disk mass on the resonance location for 
a two-planet system (without considering disk migration) 
both for the 2:1 and the 3:2 mean motion resonance. 
focused on a particular numerical difficulty 
when dealing with migrating planets. When the disk self- 
gravity is ignored, there is an inconsistency lurking in the 
torque calculation: the planet feels and reacts to the gravity 
of the gas, but the gas does not feel its own gravity. Hav- 
ing planet and disk orbit in different gravitational potentials 


causes an error in the torque calculation (see |Barureau and] 
[Masset|[2008b). [Ataiee and Kley| (2020) found that when 


this inconsistency is not corrected for, more compact plan- 
etary systems are produced. 

(2020) focused on planets 
of higher mass, capable of opening a gap. Based on the 
empirical migration speed formula developed in|Kanagawa| 
(2018b), they derive a criterion for convergent mi- 


gration of massive planets that compares favourably with 
hydrodynamic simulations. 


3.11. 


One area that was not discussed in the previous Proto- 
stars and Planets review concerns 
disk-planet interactions in binary star systems. Two basic 
configurations are possible: the planet can be in an S-type 
orbit, going round only one of the stars, or in a P-type or- 
bit, where it orbits both stars. The latter are also known as 
circumbinary planets. 

Disk-planet interactions are relevant for two aspects of 
this three-body problem. First, while the disk is present, 
the planet will interact with either the circumprimary or the 
circumbinary disk. Second, any disk in the system will be 
tidally truncated by the companion star, in a similar way 
a planet opens up a gap in the disk. A disk in an S-type 
configuration will be truncated on the outside, while a cir- 
cumbinary disk will have an inner hole. In addition, the 
binary companion can drive the gas disk to become eccen- 


tric (e.g. 1991; |Paardekooper et al.|2008), even if its 
own orbit is circular (Kley et al.|2008). 


Early work focused mostly on the circumprimary case, 


and in particular the famous case of y-Cephei (Hatzes et al. 
2003). Both Type 0 (Paardekooper et al.|2008) and Type 
(Kley and Nelson 
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VII migration 8) were investigated. 
Perhaps the most striking difference with the single star 
case is the excitation of eccentricity in the orbit of the 


planet, driven mostly by the eccentric gas disk (Kley and 
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2010). 
The discovery of Kepler 16-b by the Kepler Space Mis- 


sion (Doyle et al.|2011) has spurred research into circumbi- 
nary planets. At the time of writing, roughly a dozen cir- 
cumbinary planets are known (for a recent overview, see 
[Penzlin ef aL [2021). In all cases, the 
orbital planes of the binary and the planet are well aligned. 
One striking feature of the circumbinary planet population 
is that they orbit close to the stability limit at roughly 3.5 
binary separations (Penzlin et al.|2021). Any closer and the 
orbit of the planet eT Beene oR) At the 
same time, disk-planet interactions provide a natural stop- 
ping point for inward migration near the inner edge of the 
circumbinary cavity (Pierens and Nelson|2008). The final 
orbital parameters of the planet at the location where migra- 
tion halts depends on the physical state of the disk (viscos- 
ity, mass, thermodynamics), making circumbinary planets 
ideal test locations for disk-planet interactions (Pierens and; 


3.12. Circumplanetary Disks 


Forming planets are expected to gather material from 
the protoplanetary disks around themselves, forming their 
primordial atmospheres and potential breeding grounds of 
moons and satellites. Similar to the circumstellar disks 
around forming stars, this structure is conventionally named 
“circumplanetary disk"(CPD), although, intriguingly, it 
may not be disk-like at all. In fact, whether the CPDs are 
truly bound to their host planets is also under contention, as 
we will see later in this section. 

Simulating CPDs is numerically challenging because 
their small sizes demand resolutions much higher than typ- 
ical simulations of planet-disk interaction—orders of mag- 
nitude higher. The length scale for CPD dynamics is either 
the Hill radius ry, or the Bondi radius rp, whichever is 
smaller. For Jupiter mass planets around solar-type stars, 
ry is the smaller of the two and is about 7% of its orbit's 
semi-major axis; for Earth-like planets at 1 au, rp is stag- 
geringly small, just about 0.3% of their semi-major axes. To 
resolve a length scale so small requires special numerical 
treatment, such as grid refinement 


Wang er al. 2014 
Bailey et al.|2021), non-uniform grid geometry (Fung et al. 
Schulik et al.|2020), or local simulations (Ormel et al.| 
2015). 

Simulations allow us to extract sizes and rotation rates of 
CPDs, as well as, in models with thermal cooling included, 
the rates of gas accretion onto the planets. We will first give 
an overview of these topics, and then review the latest find- 
ings. For CPD sizes, a common misconception is that they 
equal rg/3, which may have originated from an estimate 
made by {Quillen and Trilling| (1998) through assuming the 
accreted gas has an incoming angular momentum of rz Qp. 
Numerical works have reported a range of different sizes, 
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showing that the CPD size has a more complex dependence 
on planet mass and that it is also sensitive to the thermal 
properties of the gas. Still, the value ry /3 may have some 
relevance, since it is roughly similar to the largest stable or- 


bit for ballistic particles orbiting a planet, which|Martin and| 
(2011) found to be ~ 0.4rg. 


Like its size, the kinematics of a CPD is also depen- 
dent on both planet and gas properties. Moreover, it can 
only be accurately captured in 3D simulations. Pioneering 
work done using 2D simulations have found that gas flows 
onto the CPD near the L1 and L2 Lagrange points 
[et al.[1999), but later, 3D simulations reveal that it is a ver- 


tical flow near the poles of the planets that feed the CPDs 


(D'Angelo et al.12003). In fact, in the midplane, gas gen- 
erally flows away from the planet, and the L1 and L2 La- 


grange points are exit, rather than entry, points (Machida| 
et al.]12015). The incoming angular momentum of the gas, 
which determines the final spin of the CPD, is consequently 
much different between 2D and 3D. This leads to a large 
body of work aimed at understanding the 3D rotation of 
CPDs, and some converged results are starting to emerge. 
Measuring the gas accretion rate is a main goal of CPD 
studies, since it directly informs us about the formation of 
gas giants. Simulating accretion, however, remains highly 
challenging. To do so realistically requires modeling ra- 
diative cooling, but radiative transport is computationally 
expensive. Adding to that, planets are expected to cool 
slowly, taking millions of years to cool and become gas 
giants. 3D simulations that simultaneously solve the full 
transport equation, resolve the CPDs, and track their evolu- 
tion for millions of years are simply far out of reach. Typ- 
ical CPD simulations with radiative transport employ flux- 
limited diffusion and last for tens of orbits in global models 


(Ayliffe and Bate|2009| |Szulágyi et al.|2016j|Szuldgyi|2017 
Lambrechts and Lega ; \Lambrechts et al. 
thousands of orbits in local models (Cimerman et al.|2017). 


Under the local shearing sheet approximation, the equa- 
tions governing the fluid dynamics around CPDs are greatly 
simplified. Ignoring viscosity and thermal dynamics, the 
entire set of Euler equations can be non-dimensionalized 
and be characterized by a single non-dimensional parameter 


(Korycansky and Papaloizou| 1996, |Machida et al.| 2008). 


We call this parameter the “thermal mass" of the planet: 


da = s (14) 
We distinguish it from the thermal mass of the disk, which 
equals the denominator on the right-hand side times the 
mass of the star (ie., h?M,). The thermal mass of the 
planet is then the mass of the planet divided by the ther- 
mal mass of the disk. It is straightforward to show that 
given a qth, the ratios between the key length scales ry, rp, 
and H are all uniquely determined. In particular, ry = rp 
when qth = /A/3 ~ 0.58; ry = H when qn = 3; and 
rp = H when qn = 1. Hence, planets with qu, >> 1 have 
rg > ry > H, and we call them “super-thermal” plan- 
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ets. The opposite limit is “sub-thermal” planets, which have 
qu, < land H > rg > rpg. This is a useful categorization 
that will help us identify trends in numerical results. 

Another useful categorization is whether the models 
simulate isothermal gas or not. The isothermal assumption 
was commonly employed in early work, chosen for its sim- 
plicity, but mounting evidence suggests that non-isothermal 
gas behaves vastly differently in CPDs. In the following 
sections, we will separate our review into isothermal disks 
and non-isothermal disks. 


3.12.1. Isothermal CPDs 


First, we review past findings on the rotation rate of 
isothermal CPDs. Early work by (2008); 
(2012); (2014) found rotation- 
ally supported, Keplerian disks s their planets. How- 
ever, later, simulations by L| (2015), 
and (2018) found CPDs 
POD are E aie rotating or not (nm all. Eaten ral after, 
ud Me Ee isothermal CPDs are ol eres ER E 
ter all. The main reason for the apparent disagreement is 
that they simulated planets with different values of qu, and 
made different assumptions about the physical sizes of their 
planets. Upon closer inspection, the conclusions drawn by 
different groups are in fact consistent with one another. 

Without accounting for the planet's physical size (i.e. a 
point mass), it is highly likely that CPDs are rotationally 
supported at short distances from their planets, for all val- 
ues of qu. derived empirical relations 
showing that the specific angular momentum in isother- 
mal CPDs, /cpp, scales with q when qin € 8 and q?/? 
for larger qu. In agreement with their results, 
found that [cpp ~ 0.23csrp for sub-thermal plan- 
ets. To better visualize this value, it is helpful to trans- 
late it to the “centrifugal radius" re = l2pp/(GM,), and 
in this case, we have re ~ 0.05rg. This measurement 


lcpp ^ 0.2r7,Q, for their qq, = 3 planet, which trans- 
lates to re ~ 0.013rg. These measurements are in line with 
expectations in one sense—the disk size re does appear to 
scale with min(rg, rg)—but they are surprising in that the 
disks are very small, about one to two orders of magnitude 
smaller than min(rg, rg). The left panel in Figure|9|shows 
one example of these disks. 

The fact that these disks are so small may be the reason 
why some isothermal simulations did not find rotationally 
supported CPDs. Numerical simulations of planet-disk in- 
teraction commonly include a smoothing length r, in the 
gravitational potential of the planet to avoid the singular- 
ity of a point mass. One =i interpret rs as the pe 


size of the planet. .|(2015), 
and |Kurokawa and Tanigawa E all reported ante rota- 


tion in their simulations of sub-thermal planets, and, as it 
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turns out, they all employed a rs larger than 0.05rg. Simi- 
larly, [Béthune and Rafikov\(2019) modeled the physical size 
of the planetary core as an inner boundary in local simula- 
tions, and found that the inner boundary has to be smaller 
than about 1/32 rp in order to recover the Keplerian disk. 
In reality, the physical size of the planet relative to its Bondi 
or Hill radius depends on its orbit size; it can well be a sig- 
nificant fraction of rg if the planet is close to its host star. 

Next, we review measurements of the sizes of isother- 
mal CPDs. This is a separate measurement from re, which 
is the size of the rotationally supported disk. Beyond re, the 
gas may be partially or completely supported by gas pres- 
sure, but remains bound to the planet and therefore should 
be considered a part of the CPD. On this topic, there is cur- 
rently no clear consensus on the definition of where the 
CPD ends and the circumstellar disk begins. An intuitive 
approach may be to inspect the velocity field, or, better yet, 
integrate the gas streamlines, and identify the part of the gas 
that does not move away from the planet, hence “bound” to 
it. This type of approach finds a CPD size of about i Irp 
in WME uris for sub-thermal planets 
EODEM € about 0.5rg for a i 
ets oe O eral 2008). However, this approach may not 
be ideal when coms that the gas in the CPD can "re- 
cycle" (Ormel er al 2015). 

As E gas xA flows vertically toward the 
planet near the poles and away from the planet near the 
midplane. In some isothermal simulations, this is found to 
happen everywhere 
gawa\2018), and no gas is formally bound to the planet. The 
gas near the planet is then said to be recycling, where gas 
co-orbiting with the planet flows into the CPD and back out 
within a finite time. This recycling time, which is defined as 
the enclosed mass within a given radius around the planet, 
such as the Bondi radius, divided by the outward mass flux 
at that radius, has a steep dependence on planet mass, and 
can range from one orbital period of the planet to many 
thousands of periods Fung et al.\2015). 
It is unclear how to define a CPD size in these cases. The 
concept of recycling remains important for non-isothermal 
CPDs, which we will discuss in the following section. 

Finally, we reflect on the interpretation of these isother- 
mal results. Since CPDs are generally assumed to be hot, 
dense environments, the isothermal assumption may be 
thought of as the fast cooling scenario, where heat from 
adiabatic compression is radiated away instantly, or at least 
within a time much shorter than the dynamical time of the 
gas. It is unclear whether such a limit can be reached in 
a realistic setting. Cooling is regulated by dust concentra- 
tion in the gas—too little dust and the gas cannot cool, but 
too much dust will also trap heat. Moreover, isothermal- 
ity does not only imply fast cooling, but also fast heating. 
Around an isothermal CPD, the gas is effectively receiving 
instant heating from an unknown source whenever it ex- 
pands. Such a mechanism has no clear justification, and it 
is unclear what might change if the gas cannot heat up in 
such a manner. 
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Fig. 9.— Simulations of CPDs around gin = 0.1 planets taken from|Fung et al.| (2019). The plots show slices of density with the 
planet located at the origin. Density is in base-ten logarithmic scale, in units where the unperturbed density at the planet’s location is 1. 
The isothermal simulation on the left shows a dense, flatten CPD, while the adiabatic one on the right is nearly spherically and orders of 
magnitude lower in density. Thermodynamics is critical in determining CPD structure. 


An alternative way to interpret isothermal simulations of 
CPDs is that they may resemble the end points of CPD evo- 
lution. The CPD may begin as a hot system, but after a suf- 
ficiently long time, its temperature should equilibrate with 
the surrounding circumstellar disk, assuming it has not yet 
dissipated. In this interpretation, the thermal history of the 
CPD, i.e. how it cools down over time, may play a role in 
the CPD’s final state, but is clearly not captured in isother- 
mal simulations. Despite that, some features in isothermal 
CPDs do appear consistent with observations. The giant 
planets in our solar system all host extensive satellite sys- 
tems likely born out of rotationally-supported CPDs. So 
far, rotationally-supported CPDs have only been found in 
isothermal simulations, and the r, measured are consistent 
with the sizes of the satellite systems around our gas giants 


(Fung et al.|2019). This suggests that moon and satellite 


formation may have occurred in nearly isothermal CPDs. 


3.12.2. Non-Isothermal CPDs 


Figure [9| demonstrates how sensitive the CPD structure 
is to the equation of state of the gas. Merely switching from 
an isothermal one to an adiabatic one produces a drasti- 
cally different density structure around the planet. In pur- 
suit of realistic thermodynamics, we must consider a suite 
of model parameters, including the adiabatic index, vis- 
cous heating of gas, dust-to-gas ratio, dust opacity «, and 
accretion luminosity of the planetary core Lace. In non- 
isothermal simulations, the dynamically varying disk tem- 
perature makes it difficult to define qn. Roughly, we can 
estimate h and qąn using the temperature of the disk at the 
planet's location, before the planet is introduced, given by 
h ~ ef vy where cs = \/P/p is the sound speed if the gas 
were isothermal; though, one should bear in mind that the 
temperature of the disk after the planets are introduced can 
become very different, and is dependent on radiative trans- 
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port parameters that are generally not the same between dif- 
ferent simulations. 

Simulations of non-isothermal CPDs have overwhelm- 
ingly resulted in CPDs that are supported by gas pres- 
sure with weak or no rotation at all (Ayliffe and Bate 


2009} |D’Angelo and Bodenheimer 2013} |Szulágyi et al. 
2016; Lambrechts and Lega\2017 3 


These CPDs are not disk-like, but spherically symmetric 
envelopes, making the term circumplanetary “disks” rather 
inaccurate. In some exceptions, the CPDs can become more 
disk-like. found that, for qu, = 8, if 
the temperature near the planet is fixed to a value < 2000 
K, one can recover rotationally-supported disks. |Szuldgyi 
showed that the disk morphology becomes flat- 
tened with shock surfaces forming above the and below the 
CPD when qu, Z 24. Then, in a lower mass range where 


^2 


dh = 0.1 ~ 0.5,|Lambrechts and Lega| (2017) found that 


when the product Lace exceeds a certain threshold value, 
heat transfer in the CPD becomes advection-dominated and 
the morphology of the disk becomes flattened. 

One can measure the size of these envelopes using kine- 


matics same as we do with isothermal disks. |D'Angelo 
and Bodenheimer| (2013) found that the bound envelopes 


are about 0.3 ~ 0.7rg for qu, œ 0.05 ~ 0.2, where the 
envelopes take up a smaller fraction of rg with increasing 
qth: measured ~ 0.2rp for 
their qn ~ 0.56 planets. The results of 
and(Lambrechts and Lega|(2017) appear 
in line with each other. (2017) did not iden- 
tify strictly bound envelopes, but instead defined the enve- 
lope as the region that can cool over time. They measured 
sizes of ~ 0.2rg when qu, = 0.38 — 0.75, and ~ 0.07rp 
when qu; = 1.9. The more simplified adiabatic simula- 


tions by (2019) similarly found envelope sizes 


of ~ 0.2rp for qn S 1. Overall, 0.2rg is the common 


^ 
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answer for qu, ~ 0.1 — 1 for both radiative and adiabatic 
models, which suggests that CPDs begin their evolution un- 
der roughly adiabatic conditions. For super-thermal planets, 


Lambrechts et al.\(2019) reported bound envelope sizes of 
~ O.3rg. Instead of treating radiative transport, [Kurokawa] 
(2018) modeled cooling (and heating) using 


a simple thermal relaxation treatment. They found that the 
size of the envelope depends on the thermal relaxation time. 
When the cooling time is 0.01051, it can become as large 
as rp, larger than both isothermal and adiabatic CPDs. The 
envelope size as a function of the cooling rate is therefore 
non-monotonic. 

The gas accretion rate is another important measure- 
ment one can make in radiative simulations of CPDs. This 
is a frontier of research and a unified picture has not yet 
emerged. For instance, (2019) reported 
that Jupiter-mass planets can form from 15 Earth-mass 
cores within a fraction of the disk lifetime, and showed 
that their measurements were in qualitative agreement with 


1D gas accretion models, such as those by|Piso and Youdin 
(2014) and|Lee et al.|(2014). Their models did find bound, 


pressure-supported envelopes of order the size of rp, al- 
beit smaller than commonly assumed in 1D models, which 
supports the idea that envelopes are effectively 1D. In con- 
trast, reported that 1D models 
and 3D simulations produce drastically different outcome, 
where the 3D, recycling flow of the gas prevents the planet's 
envelope from cooling. In their models, they found that re- 
cycling occurs throughout the entire envelope, all the way 
down to the core. In this scenario, none of the gas is truly 
bound to the planet, and it always has a finite time to cool 
before being recycled back to the circumstellar disk. This 
can eventually stunt the growth of the envelope and prevent 
the planet from ever developing into a gas giant. While 
and 
both simulated sub-thermal planets (q,,=0.7 and 0.18, re- 
spectively), different assumptions were made about other 
local conditions, including disk temperature and opacity. 
These factors will like be critical in understanding how 
these differing reports may be united into a single, self- 
consistent model of 3D gas accretion. 

Finally, it is worth mentioning that other physical effects, 
such as magnetic fields, can be important shaping CPDs. 
Global MHD simulations with adaptive mesh refinement 
were presented in [Gressel et al.| (2013), who found a CPD 
whose structure displays high levels of variability, and the 
generation of helical magnetic fields that launch magneto- 
centrifugally driven outflows. 


4. Tools and Techniques 


The progress that has been made in this field has been 
primarily enabled by hydrodynamic simulations. This sec- 
tion discusses briefly the basic challenges of this type of cal- 
culation 
[2018a), the state of the art, and points to the current com- 
munity codes. 
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4.1. 


Hydrodynamic simulations of planet-disk interactions 
share many challenges with hydrodynamic problems in 
other fields. One of these is the problem of a wide range 
of important scales that need to be resolved. Global disk 
evolution occurs on scales of ~ 100 AU. For a planet lo- 
cated at 1 AU, the length of an orbit, over which libration 
in the horseshoe region occurs, is 27x AU. The pressure 
scale height in the disk at this location is H ~ 0.1 AU. 
The width of the horseshoe region for a Earth-like planet is 
< H (see e.g. [Paardekooper and Papaloizou[2009b). The 
plume responsible for the thermal torque depends on the 
thermal diffusivity but is usually much smaller still 
(2017). Finally, the physical radius of an Earth-like planet 
is ~ 10-?H. It is likely that all of these scales are talk- 
ing to each other. For example, global disk evolution reg- 
ulates the inflow of solids, which in the end are responsi- 
ble for the plume governing the thermal torque. Different 
methods have different ways of trying to meet this chal- 
lenge: Adaptive Mesh Refinement for grid-based codes, 
while Smoothed Particles dynamics naturally yields high 
resolution in regions of high density. 

A second challenge arises due to the range of important 
time scales. The longest important time scale is the disk 
life time, i.e. millions of years, which is computationally 
unfeasible. Planet-disk interactions simulations need to be 
run usually until a quasi-steady state is reached. For gap- 
opening planets, this can be many thousands of orbits (e.g. 


Dempsey et al.|2021). Every orbit contains ~ 1000 time 
steps, depending on the numerical resolution 


Benítez-Llambay|2015). Each simulation then takes many 
millions of time steps, quickly making 3D computations 


very expensive. As a rough guide, a grid-based simulation, 
including only basic hydrodynamics, consisting of ~ 10? 
computational cells (i.e. ~ 1000 cells in each direction, 
see will take ~ 1000 seconds per 
time step on a modern CPU. Each simulation then takes 
~ 10° CPU hours. Additional physics in the form of radia- 
tive cooling, magnetic fields, etc, will further increase the 
computational cost. 

Since horseshoe drag relies on the material conservation 


of vortensity and entropy 199] 
2008a\ [Paardekooper and Papaloizou|2008| |2009a; 
and Masset|2009| |Masset and Casolij2009), it is important 


that numerical schemes adhere to these as much as possi- 
ble. Traditionally, staggered mesh codes, have an advan- 


tage here over those built on Riemann solvers 
Benítez-Llambay|2015). One possible explanation for this is 


that most Riemann solvers are built using directional split- 
ting and use 1D Riemann solvers, while for example vor- 
ticity is a multidimensional quantity. A truly multidimen- 
sional Riemann solver can substantially improve for exam- 
ple the lifetime of vortices in disks or MHD structures (e.g. 
[Paardekooper|[2017). SPH in principle has 
excellent conservation properties for quantities that can be 
specified on a per-particle basis (e.g. angular momentum, 
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entropy), and under certain conditions for potential vorticity 


(Frank and Reich|2003). However, since SPH simulations 


of protoplanetary disks usually result in an artificial viscos- 
ity equivalent to a ~ 107? (Price et al.|2018b), conserva- 
tion of potential vorticity in horseshoe dynamics is washed 
out by the action of viscosity. 

Another potential issue arises in particular for hydro- 
static equilibria, which could be the planet atmosphere or 
the vertical disk profile. The numerical scheme should pre- 
serve such stable states, meaning that the scheme would 
be ’well-balanced’. This is traditionally a problem for 
fractional step methods as implemented in many Riemann 
solvers, but solutions do exist (e.g. 

Shocks play an important role in particular for high- 
mass planets, capable of opening up a gap 
[and Papaloizou|1996). However even for low mass planets 
in inviscid disks shocks do appear 


[2001). It is therefore crucial that numerical methods cor- 
rectly treat shocks. Here Riemann solvers have an obvious 
advantage over other methods, since the correct treatment of 
shocks forms an integral part of the solution method. This 
allowed [Dong et al.) (2011) to obtain the correct evolution 
of a density wave excited by a low-mass planet towards a 
shock using ATHENA (Stone et al.|2008). 

As protoplanetary disk models become more sophisti- 
cated, so does the amount of physics that should be part 
of planet-disk interaction simulations. The field has moved 
from 2D and 3D, locally isothermal hydrodynamical mod- 


els [2003), to ideal MHD simu- 
lations (e.g. , to radiation- 
hydrodynamical models (e.g. 
[2008), to non-ideal MHD models (McNally et al.|2017). In 
addition, multiple phases are often included (i.e. gas and 
solids, see e.g. 
Llambay and Pessah|2018). Needless to say, the inclusion 


of all these important effects drives up the computational 
costs. 

Finally, the paradigm shift away from viscous, turbu- 
lent accretion disks towards laminar, wind-driven accretion 
disks brings a relatively new challenge. A disk where diffu- 
sion is negligible has a very long memory. This means that 
initial conditions will matter a lot more compared to pre- 
vious, viscous disk models. Planets will remember where 


they have been (see e.g.|Paardekooper|2014), which makes 


realistic initial conditions more important than ever. 


4.2. Community codes 


Below, we point to community codes that have been used 
recently to study planet-disk interactions. Some of these 
have been around since the comparison problem of. 
(2006). We mainly focus on codes that are 
freely available for download. The list is not exhaustive, 
but should give the reader an idea of the kinds of codes that 
are currently used. The feature lists for each code is incom- 
plete as well: many codes have a modular structure and new 
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modules are added on a regular basis. 


4.2.1. Staggered Mesh Codes 


We start off by describing the family of grid-based meth- 
ods that use a staggered mesh. These were termed ’upwind 


methods’ in |de Val-Borro et al.| (2006). These share their 
core hydrodynamical algorithm with ZEUS (Stone and Nor-| 


1992) and are based on finite differences. 
FARGOG3D: The original version of FARGO 


was specifically aimed at planet-disk interactions, 
and its successor FARGO3D 
has inherited the features that make it straightforward 
to set up a planet-disk interaction problem. FARGO3D can 
do MHD and can run on GPUs. It has a module to in- 


clude multiple dust species (Benítez-Llambay et al.12019). 
Recent applications to planet disk interactions include 
ally et al.|(2020a) and |Debras et al.| (2021). 


Other well-known staggered mesh codes include the 


original NIRVANA code (Ziegler and Yorke| 1997), which 


was never publicly released, but was part of the 2006 com- 


parison problem in several incarnations 
[2006), the FARGOCA version of FARGO maintained in 
Nice (Lega et al.[2014), and the code used in|Vorobyov and| 
OTS. 


4.2.2. Riemann solver Codes 


Many of the current community codes have Riemann 
solvers at their core. We list them here in alphabetical order. 
ATHENA ++: This is the successor of the Athena MHD 
code (Stone et al.|2008), ported to C++ and with new fea- 
tures added (Stone et al.|[2020). It has a Riemann solver 


at its core, and is fourth order accurate in space and time 


in regions of smooth flow (Felker and Stone|2018). The 
code was first used to simulate planet-disk interactions by 
(2015b) and (2015). More recent 
publications include (2021) and[Kuwahara and] 
[Kurokawa] (2020) 

NIRVANA: The public version of NIRVANA uses Rie- 


mann solvers (Ziegler|2004), and has support for non-ideal 


MHD, adaptive mesh refinement and self-gravity. It has 


been used recently in McNally et al.|(2017). 
PLUTO: The finite volume code PLUTO 


uses a Riemann solver at its core, and has 
support for MHD and solid particles. A version that runs 
on GPUs was presented in{Thun et al.|(2017). Examples of 
recent work on planet disk interactions include 
and [Penzi eral 2021) 

Other, not publicly available codes using Riemann 
solvers include RODEO 
and PEnGUIn (Fung et al.|2014). 


4.2.3. High-order Finite Difference Codes 


This class of methods achieves high accuracy by taking 
high-order finite differences. The resulting scheme needs to 
be stabilized by artificial viscous terms. 
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PENCIL: This is a high-order finite difference MHD 


code (Pencil Code Collaboration et al.12021) and has sup- 
port for solid particles in addition to gas (Youdin and Jo- 


2007). For recent applications on planet-disk inter- 
action see (2016) and|Yang and Zhu|(2020). 


4.2.4. Smoothed Particle Hydrodynamics 


SPH methods are popular for disk simulations, since 
they naturally allow for complex geometries such as warps, 
and because of their oct-tree data structure, the inclusion of 
self-gravity is straightforward. 

PHANTOM: This Smoothed Particle Hydrodynamics 
(SPH) code (Price et al.|2018b) has support for MHD and 
dust e erar m) and is widely used to 
simulate accretion disks. Recent publications on planet disk 
interactions include and 
(2020). 

GADGET: The publicly available version GADGET- 
2 has support for self-gravity. For a re- 


cent application of GADGET to disk-planet interaction see 


Humphries and Nayakshin\(2019). 
SEREN: This SPH code was developed for star and 


planet formation simulations (Hubber et aL.|2011a]b). It 
supports self-gravity and cooling through radiation. A re- 
cent application to disk-planet interactions is 
lend Tsu COIS) 


For other, not publicly available SPH implementations 


see for example Fletcher et al. (2019). 


4.2.5. Other methods 


DISCO: The moving-mesh MHD code DISCO 
employs a moving-mesh approach utilizing a dy- 
namic cylindrical mesh that can shear azimuthally to fol- 
low the orbital motion of the gas, in combination with a 
Riemann solver. It is specifically designed to deal with the 
shearing motion present in astrophysical disks. Recent pub- 
lications on planet-disk interactions include [Duffell] (2020) 


and |Duffell and Chiang (2015). 
High-order discontinuous Galerkin (DG) methods (e.g. 


are slowly finding their way 
into astrophysical fluid dynamics. An implementation for 
planet-disk interactions was presented in [Velasco Romero| 
(2018). While the use of DG methods in planet- 


disk interactions is still in the exploratory phase, we have 
included them here as because of their combination of high- 
order accuracy with shock-capturing capabilities, they may 
well represent the future of this type of calculations. 
Finally, it is worth mentioning the moving mesh code 


AREPO (Springel|2010), which has been used on the disk- 


planet problem in (2014), and the meshless 
code GIZMO (Hopkins|2015), which was part of the com- 


parison paper|Fletcher et al. 


5. Modeling Observational Signatures 


In this section we introduce recent work on modeling ob- 
servational signatures of disk-planet interactions. Advance- 
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ments in theory and numerical simulations in the last decade 
play a crucial part in the emergence and development of this 
new field. More detailed discussions can be found in the 
chapter led by Bae et al. We focus on dust observations 
(ONIR scattered light and mm/cm continuum emission). 
Planet-induced gas kinematics signatures are discussed in 
more details in the chapter led by Pinte et al. 
Observational signatures of planet-induced structures 


are generally modeled in a two-step process (e.g., 
2015a). First, hydrodynamics simulations are carried 


out to calculate the spatial distribution and kinematics of 
gas perturbed by one or more planets. Dust particles of 
various sizes may be added into the simulations to obtain 
the spatial distribution of dust as they dynamically interact 
with the gas. In the second step, the resulting hydro disk 
models are fed into radiative transfer simulators, which pro- 
duce synthetic images of various kinds. The hydro models 
may be intrinsically 3D, which provide volume density and 
velocity distributions, or 2D, in which case they are often 
“puffed up" in the vertical direction to obtain 3D structures. 
The resulting synthetic images may be post-processed to 
achieve realistic angular resolutions and sensitivities, to 
facilitate direct and quantitative comparisons with observa- 
tions. 


5.1. Generic modeling of planet-induced structures in 


observations 


Planets on stationary, coplanar, and circular orbits are 
expected to produce disk structures, including spiral arms, 
gaps, vortices, and circumplanetary disks. Inclination or 
eccentricity in the planets' orbits and planet migration can 
affect the appearance of these features, or even produce ad- 
ditional ones. 


5.1.1. Spirals 


One planet can drive multiple spiral arms both inside and 
outside its orbit, as different modes in the resonances ex- 


cited by the planet constructively interfere (Bae and Zhu 
2018aj |Miranda and Rafikov\\2019a). The number and 


strength of the spirals depend on the planet mass (Mj) as 
well as disk properties such as aspect ratio (h/r), and ther- 
mal relaxation rate 
(2020). At tens of AU, spiral arms excited by planets more 
massive than Saturn may be traceable in ONIR scattered 
light using current generation of instruments (Dong and; 
2017a) at a variety of viewing angles 
2016b). As the planet mass increases, the strength of spi- 
rals in scattered light and the separation between the pri- 


mary and secondary spirals increase (Fung and Dong\2015 
Dong and Fung|2017a;|Bae and Zhu|2018b); such quan- 


titative correlations can be exploited to dynamically con- 
strain the planet masses. Specifically, a planet with My ~a 
few x M, (h/r)? (usually supra-Jovian) can drive a pair of 
prominent arms in a near symmetric configuration inside its 


orbit (Dong et al. [2015b). 


Spirals may also be detected as temperature perturba- 
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tions in disks with non-isothermal EOS as gas been com- 


pressed and decompressed crossing spirals (Bae et al.|2021 
Muley et al.|2021). Inside spiral arms, significant vertical 
motion is present and may be traceable (Bae et al.|2021). 


Because of this, modeling observational signatures of spi- 
rals generally requires 3D disk-planet interaction simula- 
tions; puffing up 2D hydro models assuming hydrostatic 
equilibrium would fail to capture the vertical motion and 
misrepresent the gas density structure, producing scattered 


light signatures that are too weak (Dong and Fung|2017a). 
5.1.2. Gaps 


Planet-opened gaps appear as annular structures with 
low intensities in observations. Under a modest disk vis- 
cosity (a = 107°), usually only one gap at planet's orbit 
forms. The radial profile of gas surface density, as well as 
the depth and width, of such a gap may be parametrized us- 
ing Mp, h/r, and a (Fung et al.|2014| 
2017). In ONIR scattered 


light, the contrast of a single gap with the adjacent disk 
can be analytically expressed using these parameters and 
the angular resolution of the observations 
[2017b). Gaps are also detectable in mm/cm dust continuum 
observations that probe the spatial distribution of mm/cm- 
sized dust. In those observations, the adjacent dust rings 
may become broad due to dust feedback 
[2018a), and have finite vertical thickness due to planet- 
driven meridional turbulence (Bi et al.|2021). Aerodynami- 
cal coupling between the gas and dust leads to dust concen- 
trations at gas pressure peaks, resulting in higher contrasts 
in dust emission observations than in scattered light 
et al.|2016). Models of the observational signatures of gaps 


are generally insensitive to 3D dynamics (Fung and Chiang 
2016} |Dong and Fung|2017b), but 3D simulations are still 


required to capture the meridional motion. 

In theory, the dissipation of each spiral arm excited by a 
planet can lead to the opening of an individual gap, thus 
one planet may open a few gaps as it drives a few spi- 
rals [Bae et al.[2017). Such signatures 
are most prominent when there is no mechanisms to dissi- 
pate spiral arms other than their intrinsic nonlinearity, e.g., 
in disks with low viscosity (a < 1074; 
[Rafikov|2001). The density waves and the separation be- 


tween these gaps are regulated by Mp, h/r, and thermal 


dition, if the planet is migrating, the spacing and morphol- 
ogy of the gaps can be further modified (Meru et al.|2019 
[Nazari ei al [2019} [Pérez er al [2019 2015}, 


and more gaps may be produced if the planet migrates in a 


non-steady way (Wafflard- Fernandez and Baruteau|2020). 


Finally, a low mass planet may open gaps in dust of cer- 
tain sizes even when it does not open prominent gas gaps 


(Dipierro and Laibe|2017). 
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5.1.3. Vortices 


Vortices can be excited by the Rossby wave instability at 
the edge of planet-opened gaps. They trap dust particles and 
can appear as emission clumps in thermal continuum ob- 
servations (Zhu and Stone{2014). The location of the dust 
concentration depends on the dust size, thus observations 
at different wavelengths may detect azimuthal variations 


in the vortex location 
land Zhu[2016). In addition, vortices produce spiral arms, 
though they are typically too weak to be detectable 
erat poro). 

When the dust-to-gas mass ratio reaches order unity at 
the gap edges, dust feedback enhances the sharpness of the 
edges, and leads to instabilities and the formation of tiny 
vortices with highly concentrated dust in continuum obser- 


vations (Huang et al.|2020; Hsieh and Lin|2020). 


5.1.4. Signatures of planets on eccentric and/or inclined 
orbits 


A massive planet on an inclined orbit may break the in- 
ner disk out from the outer disk, and cause the two parts 
to precess at different timescales (Zhuj2019). Consequently 
the two parts of the disk may become warped or misaligned. 
The inner disk may cast shadows on the outer disk, and the 
radially varying disk inclination may be visible in gas kine- 
matic observations, in a way similar to the effect of a mis- 


aligned binary (Facchini et al.|2018). A giant planet on an 
eccentric orbit can carve out a gap as wide as its orbit, much 


wider than its circular counterpart (Muley et al.|2019). Such 


a gap can also appear much broader in C/^O than in !?CO 


observations (Baruteau et al.|2021). 


5.1.5. Signatures of planets in disk kinematics 


Planet-induced kinematic signatures have been modeled 
using hydrodynamics simulations and synthetic gas obser- 
vations. showed that the rotation of a 
CPD around a giant planet produces a unique signature in 
disk channel maps well separated from the circumstellar 


disk. (2018) showed that a planet may create 


localized “kinks” in the isovelocity maps of its host disk. 


(2018) quantified the radial modulation in the 


circumstellar disk rotation due to the pressure gradients at 


gap edges. (2018) investigated gas kinematics 


caused by pressure gradients at gaps, spiral wakes, and vor- 
tices. and showed 
that a planet opening a deep gas gap can drive transonic ver- 
tical motions inside the gap potentially detectable in disk 
Moment 2 maps. More discussions on these topics can be 
found in the chapter led by Pinte et al. and in the recent re- 


view article by Disk Dynamics Collaboration et al.| (2020). 


5.2. Planet-disk interaction as the origin of observed 
disk structures 


Hundreds of protoplanetary disks have been imaged by 
the latest generation of instruments. In most of them, struc- 
tures that resemble planet-induced have been found. Obser- 
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vations are now being quantitatively compared with disk- 
planet interaction models to infer the presence of planets 
forming in disks, and to constrain their properties, such as 


mass and orbit (Bae et al.|2018| Zhang et al.|2018 
[et _al.|/2019). These planets are extremely difficult to find 
directly (e.g.,|Keppler et al.|2018); yet they are the key- 


stones in testing planet formation theories, as they provide 
direct constraints on the location, timescale, and local envi- 
ronment of planet formation. 

There have been numerous models where disk features 
are explained with planet-disk interactions. Good quantita- 
tive agreement between theory and observation in many of 
them supports the planetary origin of those features. Below 
we review some representative examples focusing on struc- 
tures seen in images. Disk-planet interaction models have 
also been tested using other observational diagnostics, e.g., 
the spectral energy distribution of the disk (e.g., in the case 


of CI Tau, Muley and Dong|2021). 


1. HL Tau. As the first protoplanetary disk targeted by 
ALMA with its long baseline configuration 
[Partnership et al.|2015). the HL Tau disk attracted 
immediate attention when several gaps at tens of AU 
were revealed. Several groups attempted to explain 
its gaps using planets. Assuming a modest viscosity 


of a ~ 107%, (2015a), 
(2015) and Jin er al.| (2016) reached the same con- 


clusion that the three main gaps can be explained as 
opened by three Saturn mass planets at around 13, 33, 
and 70 AU (Figure[I0]. 


HL Tau (ALMA Partnership et al. 2015) 


Observation 


Simulations 


Jin et al. 2016 


Dipierro et al. 2015 


Dong et al. 2015 


Fig. 10.— Top: ALMA observations of the HL Tau disk, show- 
ing multiple rings and gaps in dust continuum emission at mm 


wavelengths (ALMA Partnership et al.12015). Bottom: Simula- 


tions by three independent teams that show the three main gaps in 


the disk may be explained as opened by three Saturn mass planets, 

if the disk has a modest viscosity o, ~ 10? (Dong et al. 2015a 
2015 2016). See section|5.2|for details. 

2. MWC 758. The MWC 758 disk was one of the first to 

be discovered to harbor a pair of near-symmetric spi- 


ral arms in scattered light (Grady et al.|2013 
2015). COTSb) and 
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showed that a pair of such arms can be ex- 
cited by a super-thermal mass planet on the outside 
at ~100 AU. Later, ALMA mm continuum observa- 
tions showed that the disk harbors a cavity ~ 40 AU 
in radius, and in the outer disk there are two emis- 


sion clumps at ~ 50 and 80 AU (Boehler et al.|2018 
2018c). The northern clump has a coun- 


terpart in cm emission revealed by the VLA obser- 
vations, but not the southern clump 
[2019). explained these addi- 
tional features by involving a second planet ~ 2M; 
in mass at 35 AU. Together with the outer planet, they 
each open a deep gap, and trigger the formation of a 
vortex trapping mm-sized dust at the gap edge (inner 
edge for the outer gap and outer edge for the inner 
gap). The southern clump is in the process of decay- 
ing, thus it does not trap cm-sized dust well. MWC 
758 is an excellent example of using disk-planet in- 
teraction models to explain disk structures in multi- 
wavelength observations (Figure[11). 


MWC 758 


1.0 um (VLT/SPHERE) 1 mm (ALMA) 1 cm (JVLA) 
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Fig. 11.— Top row: Three observations of the MWC 758 disk 
at different wavelengths on the same scale. From left to right: a 
scattered light image at ~ 1um ,an ALMA 
continuum emission map at ~ 1mm (Dong et al.|2018c), and a 
JVLA continuum emission map at ~ 1cm (Casassus et al.|2019). 
Bottom row: Three synthetic observations of the same kind as the 
panels above (Dong et al. D015b][Baruteau et al.[2019). The orbits 
of the planets in the models are traced out by the green dashed 
circles. There is one planet in the model on the left, and there are 
two planets in the model shown in the middle and right panels. 
The models match the observations reasonably well. See section 


[5.2]for details. 


3. HD 169142. (2019) discovered a double 


gap and triple ring structure in a compact configu- 
ration at around 70 AU in the HD 169142 disk us- 
ing ALMA mm continuum observations, and showed 
that the structure can be produced by an inwardly mi- 
grating ~10 Earth mass mini-Neptune at ~ 68 AU in 
a disk with low viscosity. This is one of first cases 
in which migration is required in order to explain ob- 
served disk structures, and the predicted planet mass 
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is among the lowest (Figure [12}. 


A6 | arcsec 
o 
o 


0.0 0.0 
Aa / arcsec Aa | arcsec 


Fig. 12.— Left: ALMA 1.3 mm image of the HD 169142 disk. 
Right: Synthetic observations of a disk with a migrating mini- 


Neptune planet of 10 Mg (Pérez et al.|2019). See section[5-2]for 


details. 


While in the above examples, disk-planet interaction 
models reproduced real observations, considerable uncer- 
tainties still exist. First of all, the properties and even the 
number of planets needed to explain certain structures de- 
pend on disk properties, some of which remain poorly con- 
strained. For example, in the case of HL Tau [Bae et al.] 
showed that most observed gaps can be explained 
by a single ^ Neptune mass planet if the disk viscosity if 
lower than a = 1074. A similar example has been shown 
in the case of AS 209, in which 5 gaps may be produced 
by one planet (Zhang et al.[2018). Further discussions on 
this topic can be found in the chapter led by Bae et al.. Sec- 
ondly, alternative theories that do not involve planets may 
also reproduce observed structures. Whether these mecha- 
nisms are operating in a significant fraction of disks is still 
under debate. For example, in the case of MWC 758, the 
gravitational instability (GI and 
central binaries (with a stellar or brown dwarf secondary) 
have both been proposed as alternatives to planets in re- 
producing some of the observations — 
shows that GI can reproduce the observed symmetric pair of 
spiral arms (see a similar example in the case of Elias 2-27; 


Meru et al.|2017), and (2020) used a central 


binary to explain the spirals as well as a few other features 
(sce also [Price etal [018a] Pablere er ar [2019] 2020). 
The most direct way to verify disk-planet interaction as 
the origin of observed disk structures is to find the planets 
predicted by the models. The best example here is PDS 
70, in which a wide gap ^ 60 AU in radius has been 
found in scattered light and dust continuum observations 


(Hashimoto et al. [2012] [Dong et al. [2012] [Hashimoto et al. 
2015). (2011/|2012b) showed that multiple plan- 


ets may open a giant common gap much wider than the in- 
dividual gap opened by a single planet. Together with dust 
filtration at the outer gap edge, such giant cavities, a.k.a. 


transitional disks, can be explained (see also |Duffell and 
[Dong|2015] [Dong and Dawson|2016). Two giant planets 
have since been found inside PDS 70's cavity (Keppler et al. 
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2018} [Wagner er al. 20183] [Hafferr et al [2019 
2020), and models show that they can well reproduce 


the observed disk structure (Bae et al.|2019). 

However, the success of PDS 70 is a rare case at the 
moment. In general, testing the planetary hypothesis for 
observed disk structures has been challenging. Planet can- 
didates have been found in a few other systems and pro- 
posed to be the driver of observed disk structures, such 


as MWC 758 (Reggiani et al.|2018 2019), 
LkCa 15 (Sallum et al.|2015), and HD 100546 
et al./2014} [Currie et al.|2015). However their status re- 


main uncertain at the moment. Direct imaging surveys have 
been carried out to systematically search for planets in disks 


(e.g., 2021). They remain largely in- 


conclusive regarding the presence of the feature-producing 


planets (e.g., 2018b), partially due to the de- 


generacy in the planetary luminosity between hot and cold 
start models (Fortney et al.|[2008). The difficulty in find- 
ing these predicted giant planets may be expected, as they 
may experience episodic accretion and are only detectable 
in a small fraction of their formation lifetime (Brittain et al. 
(2020). The formation of some of these predicted planets, 
e.g., the multiple Saturn mass planets in the HL Tau disk 
at tens of AU and the supra-Jovian planet in the MWC 758 
disk at over 100 AU, may also challenge our theoretical un- 
derstanding on how planets form. The comparison between 
planet demographics at birth with those around billion years 
old stars obtained from radial velocity and transit surveys 
may also offer clues on the orbital evolution of planets. 
Due to the difficulties in directly confirming the theory- 
predicted planets in disks, indirect approaches have been 
explored. For example, a brown dwarf or a stellar mass 
secondary can drive a pair of prominent arms in a sym- 
metric configuration inside its orbit, in a way similar to a 
giant planet. In the case HD 100453, such an outer stellar 
companion (with a primary/secondary mass ratio of roughly 
10:1) and a pair of spiral arms have been found (Chen et al.| 
2006} [Wagner et al. 2015} 
2019). The companion’s orbit has been traced out 
(Wagner et al.|2018b), and models and analysis have shown 
that the companion is likely responsible for driving the spi- 
rals 
[2020). The HD 100453 case is direct proof of an external 
companion exciting a pair of spirals as predicted by theory. 
The pattern speed in asymmetric structures may also be 
used to test their disk-planet interaction origin. This has 
been attempted in the cases of MWC 758 and SAO 206462, 
both of which show a pair of symmetric arms (Benisty et al.| 
2015 
2016). Theory predicts that if the spirals are produced by a 
companion, the spirals and the companion should orbit the 
star at the same angular velocity (assuming the companion 
is on a circular orbit). measured 
the pattern speed of the two spirals in the MWC 758 disk, 
and concluded that they are consistent with being excited 


by an external companion at ~ 170 AU. (2021) 
applied the same technique and measured the pattern speed 


Paardekooper, Dong, Duffell, Fung, Masset, Ogilvie, Tanaka 


in the SAO 206462 disk. They found that if the two spirals 
are co-moving, they are consistent with being excited by a 
planet on the outside at ~ 86 AU. Meanwhile, current data 
does not provide strong constraints on whether they are co- 
moving or not, and if they are independently moving, they 
may be excited by two planets at ~ 120 and 50 AU. 


6. Summary and Conclusions 


Despite having been around for more than 40 years 
|dreich and Tremaine|1979), the theory of disk-planet inter- 
actions is still seeing significant progress. This is driven in 
part by an ever improving of our understanding of the intri- 
cate processes leading to angular momentum exchange be- 
tween planet and disk (for example thermal and dynamical 
torques, and the impact of solids) and in part by a better un- 
derstanding of protoplanetary disks, where a paradigm shift 
has occurred from viscous to wind-driven accretion disks. 

Since the last edition of Protostars and Planets, planet- 
disk interactions has become a subject that can be observed 
directly, rather than inferred from the current exoplanet 
population. We can use the theory of planet-disk interac- 
tions to interpret observations (see section 5}. and a partic- 
ular exciting prospect for the coming years is observation- 
driven improvement of the theory (e.g.|Nazari et al.|2019). 
This is necessary, since despite all recent progress, the prob- 
lem of migration has not been ’solved’. If history is a 
guide, introducing new physics into the problem produces 
different torques that need to be understood. On the nu- 
merical side, given the small length scales involved in for 
example thermal torques, clever use of adaptive resolution 
may provide a way forward to understand how these small 
scales interact with the largest scales that in the end will 
determine how much mass can potentially flow towards the 
planet. The same holds for circumplanetary disks: the ex- 


citing observation of |Benisty et al.| (2021) means we need 


the most realistic CPD models possible, and probably re- 
solve the planet at the same time (Zhu et al.[2021). 

In view of all this, and the current influx of detailed 
observational data on protoplanetary disks and embedded 
planets, we expect again a lot of progress in our understand- 
ing of planet-disk interactions over the next few years. 
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